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PREFACE 


This document, prepared by the Dynamics and Loads Section, Martin Marietta Corporation, 
Denver Division, under Contract NASS- 1 1996, presents the results of a study for the purpose 
of developing a computer program system for dynamic simulation and stability analysis of 
passive and actively controlled spacecraft. The study was performed from May 1973 to 
April 1975 and was administered by the Goddard Space Flight Center, National Aeronautics 
and Space Administration, Greenbelt, Maryland, under the direction of Joseph P. Young. 

Upon delivery, the computer program and associated documentation was checked in detail 
by Harold P. Frisch. This document incorporates both the original Martin work and the 
supplementary material prepared at the Goddard Space Flight Center. 

The digital computer program, DISCOS (Dynamic Interaction Simulation of Controls and 
Structure), has been extensively annotated and tested on a range of problems that should 
have exposed nearly all theoretical errors and programming bugs. 

From its inception in 1973, DISCOS has been designed to grow as new needs and more 
efficient computational techniques develop. This feature makes it impossible to define a 
final version. To circumvent this problem, the official release version will contain only 
those additions to the delivered program that enhance program documentation and user 
interface capability and correct programming errors. 

Included in this version are more than 10,000 comment cards and a capabihty to routinely 
direct the computer to output on the line printer virtually all computation along with 
explanatory alphanumeric statements. A large percentage of the comment cards are in 
subroutines DEFl, DEF2, . . ., and DEF5. These subroutines are composed entirely of 
comment cards and provide the user with an area in the source file for keeping documenta- 
tion current. In particular, subroutine DEF5 contains a narrative description of the program 
and its current capabilities. 

For the uninitiated reader, it probably would be best to speed-read subroutine DEF5 to 
obtain a quick overview of the capabilities of the program and the solution techniques 
applied before reading this document. 

The potential user should be aware of the fact that DISCOS is not intended for simple prob- 
lems. It is primarily for problems which were heretofore intractable. Consequently, the 
theoretical basis for the program is highly advanced, and the computation algorithms are 
designed for the efficient processing of the equations associated with large multidegree-of- 
freedom systems. 
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As an aid to the user, the paper on which the derivation of the coupled .flexible-body equations 
of motion is based and a paper that outlines the solution method and comments on results 
obtained from several DISCOS applications appear as reference material at the end of Volume 
I. Volume I contains all relevant theoretical work. 

Volume II describes the DISCOS program and its support programs. The user is encouraged 
to refer to the comment cards provided in all DISCOS subroutines for additional descriptive 
information. The comment cards found in the DISCOS subroutines are intended to provide 
a link between the computer code and the theoretical equations provided in Volume I. 

To effectively interface with the program, the user must be able to write the subroutines 
that will define all non-gyroscope forces and torques. The user is provided with a clean 
interface. When the load vector associated with the effects of springs, dampers, motors, gas 
jets, etc. is defined and properly stored in the computer memory, DISCOS will automatically 
transform it into the appropriate generalized form required by the formulation. 

The inclusion of effects such as aerodynamic loading and thermodynamic deformation is 
more difficult. However, the methodology is analogous to that used for including gravity- 
gradient effects. 

The methodology for including loads associated with springs, dampers, motors, gas jets, 
constraints, etc. is found in Volume II and in the comment cards of the appropriate sub- 
routines referred to in Volume II. 

DISCOS is probably the most powerful computational tool to date for the computer simu- 
lation of actively controlled coupled multiflexible-body systems. It is not an easy program 
to understand and effectively apply, but it is not intended for simple problems. The user 
is expected to have an extensive working knowledge of rigid- and flexible-body dynamics, 
flnite-element techniques, numerical methods, and frequency-domain analysis. 
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A DIGITAL COMPUTER PROGRAM 
FOR THE DYNAMIC INTERACTION SIMULATION OF 
CONTROLS AND STRUCTURE (DISCOS) 
VOLUME I 

Carl S. Bodley, A Darrell Devers, and A. Colton Park 
Martin Marietta Corporation 
Denver. Colorado 

and 

Harold P. Frisch 
Goddard Space Flight Center 
Greenbelt, Maryland 


I. INTRODUCTION 

Modem society has derived numerous and varied benefits from Earth-orbiting satellites. 

These benefits include prediction of extreme changes in weather, exploration of mineral 
and water resources, materials research, medical research and experimentation, solar-system 
studies, and increased communications and defense capabihties. The design and ultimate 
development of these satellites requires extensive analytical and experimental studies to 
ensure complete confidence in the overall ability of the total system to perform its required 
functions. 

Two of the most important and potentially most difficult of these studies are the analysis 
of dynamic response and the prediction of stability characteristics. In recent years, the 
National Aeronautics and Space Administration (NASA) and members of the aerospace 
industry have expended much effort in analyzing these phenomena. Although these are 
worthy efforts, most have been somewhat limited in scope because of the extreme com- 
plexity of the problems. Factors that make these studies so complex are: (1) the imposition 
of a spin rate for either stabilization or artificial gravity ; (2) the dynamic interaction of this 
spin rate with large flexible appendages; (3) the complex environmental loadings, including 
gravity gradient and solar radiation; and (4) the intricate control logic required for. maintain- 
ing stability or for executing orbital maneuvers. 

The authors at Martin Marietta Corporation acknowledge the assistance provided by Goddard 
Space Flight Center (GSFC) personnel. Harold Frisch and James Donohue contributed 
valuable technical comments and suggestions throughout the program. In particular, they 


developed the basic approach to be used for data input, defined exactly what should be in- 
cluded in the transfer function and stability analysis portion of the program, and defined 
eight of the eleven demonstration problems. They provided the associated data that were 
used to verify both the nonlinear time response and linear transfer function and the stability 
analysis portions of the program. Raymond Welch provided the subroutine used to generate 
root-locus plots. Dr. William Case provided valuable advice on interfacing with NASTRAN 
output and generated the demonstration problem used to validate the interface subroutine 
(NASFOR).* During the early program development stage. Dr. James Mason offered signifi- 
cant advice on the need to compute internal forces at the interconnect points. Reginald 
Mitchell contributed invaluable advice on requirements for making the program compatible 
with the GSFC IBM 360/95 computer system. In addition, he supplied the contractor with 
a 360-compatible plot package, furnished the contractor with a self-authored subroutine for 
reading NASTRAN output, and was responsible for running all demonstration problems on 
the 360/95 computer. Finally, the authors acknowledge the encouragement and efforts of 
Joseph P. Young, Technical Monitor, who made numerous valuable comments and suggestions 
throughout the study. 

A. Overview 

The state-of-the-art dynamic response analysis of a system of connected bodies is currently 
restricted to topological systems of connected rigid bodies with (possibly) flexible terminal 
bodies. Because of the complex orbiting configuration and mechanical systems proposed 
for future space programs, the limitations of the current technology are severely restrictive. 
This document presents a more comprehensive formulation that is capable of considering 
any body of the total system as flexible and that is not restricted to a specific connection 
arrangement. 

Applications of such methods and program systems are numerous and include simulation of 
the Space Shuttle payload deployment/ retrieval mechanism, solar-panel-array deployment, 
antenna deployment, analysis of multispin satellites, and analysis of large, highly flexible 
satelhtes. 

This approach provides a general-purpose modeling capability for dynamic simulation and 
stability analysis of passive and actively controlled spacecraft. In particular, the following 
items are considered: 

• Time-domain solution of the nonlinear differential equations of motion that 
describe total or nominal response^ of the complete spacecraft system idealized 
as a collection of interconnected flexible (or rigid) bodies 

• Linearization of the governing differential equations by numerical means 


♦The NASFOR subroutine has been superseded by a special NASTRAN DMAP program and associated preprocessor 
program written by Harold P. Frisch. 

fin certain cases, the total response of the dynamic system may be considered to be equilibrium state motion (nominal 
response) plus perturbation motion with respect to the equilibrium state. 
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• Time-domain solution of the linearized differential equations that describe the 
perturbation response of the complete spacecraft system about some predetermined 
(calculated or user-specified) nominal motion 

• General frequency-domain stability analysis corresponding to the linearized space- 
craft representation 

• Provision for arbitrary (explicitly time-dependent) loadings and environmental 
interaction, such. as gravity gradient and thermally induced deformations resulting 
from solar radiation 

B. Description of the Physical System 

The physical system undergoing analysis may be generally described as a cluster of contiguous, 
flexible structures (bodies)' that comprise a mechanical system such as a spacecraft. The 
entire system (spacecraft) or portions thereof may be either spinning or nonspinning. Member 
bodies of the spacecraft are capable of undergoing large relative excursions such as those of 
appendage deployment or rotor/stator motions. The general system of bodies is, by its 
inherent nature, a feedback system in which inertial forces (such as those due to centrifugal 
and Coriolis acceleration) and the restoring and damping forces' are motion-dependent. Also, 
the system may possess a control system in which certain position and rate errors are actively 
controlled through the use of reaction control jets, servomotors, or momentum wheels. 

Bodies of the system may be interconnected by linear or nonlinear springs and dampers, by 
a mechanism that is a combination of gimbal and slider block, or by any combination of these. 
Also, any two bodies of the system may be free (one from the other) and possess six degrees 
of relative motion freedom. Several or all of the six degrees of relative motion freedom be- 
tween two bodies may also be a prescribed function of time (including zero motion). 

For further introduction of the physical system, consider an Dlustrative example, such as 
the system of bodies of figure 1, and let this example be the prototype for all subsequent 
discussion and development. 

In figure 1, a nontppological tree configuration has been deliberately indicated. There are 
five hinges and four bodies, thus one closed path. Consecutive integer labels ate used for 
body reference points, hinges, sensor points,* and inomentum wheels. The numerical 
order iii each of the four label sets may be random; however, it is understood that body 1 
(which may be any body of the system) is sasociated with hinge 1. ■ 

For each body of the system, there is a body-fixed, right-handed reference frame, whose 
origin may be at the body’s mass center or at some structural hard point on the body. (A 
body’s elastic deformation is me^ured in its reference frame.) 


*A “sensor poinf ’ is any point at which kinematic data must be obtained (e.g., where an attitude sensor is located). 
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Figure 1. Labeling scheme for example system. 

In this document, a hinge is defined as a pair of structural hard points (see figure 2) with a 
point situated on each of two contiguous bodies. In figure 2, point p and point q comprise 
a hinge. Clearly, a typical body may contain one or more hinge points, but a hinge may be 
associated with only two bodies. Hinge 1 is given special consideration. It is also a pair of 
points; but one of the pair is coincident with the reference point of body 1, and the other 
point of the pair is coincident with the inertial origin. Thus, motion “across the hinges” 
is used to represent system motion. The reference point of body 1 is located with respect 
to the inertial origin by an inertially referenced position vector. The attitude of the refer- 
ence frame of body 1 , with respect to the inertial frame, is represented by three Euler 
angles. Thus, six position/ attitude coordinates are associated with hinge 1. 

Each of the remaining hinges is considered in a manner somewhat similar to that of hinge 1 , 
Referring to figure 2, note that there is an orthogonal reference frame attached to point p 
and another to point q. The triad of point p may have a “natural” (or undeformed) mis- 
alignment with respect to the triad of body point m, plus additionalmisalignment caused 
by elastic deformation. The same relationship is tme concerning points n and q. 

Now six relative position/attitude coordinates are associated with the hinge of points p and 
q. Point q is located from point p with a p-frame referenced position vector. The attitude 
of the q-frame with respect to the p-frame is represented by three Euler rotations. Thus, 
if NH is the number of system hinges, 6 X NH position coordinates are to be used in 
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Figure 2. Typical contiguous bodies of the system. 

conjunction with modal displacement coordinates for defining the system’s position state. 
Note that only the time- variable position coordinates of the 6 X NH set of candidates are 
considered as state-vector elements. (The position coordinates whose rates are constrained 
to zero are not included; however, the position coordinates themselves need not be zero.) 

The system of bodies usually has a number of so-called sensor points. A sensor point is 
defined as a structural hard point with an attached right-handed orthogonal reference frame 
that is used for a variety of purposes. A sensor point may be used to enable the program 
system to monitor the position, attitude, or the rates associated with a specific structural 
hard point. For example, a rate gyro, a star tracking device, or another motion/position 
sensing device is physically situated at a sensor point. Also, a sensor point is used as a point 
of application of a force or torque vector (see figure 2). 

The system of bodies may contain buUt-in momentum wheels, of which some are constant 
speed wheels and others are variable speed. The variable speed momentum wheels are motor 
driven; the shaft torque results from a given control law. Each momentum wheel of the 
system must be associated with a sensor point because, for a general flexible body, the gyro- 
scopic coupling is influenced by elastic motion. 

As is indicated in figure 1, the system may be in a nontopological tree configuration. The 
methods employed in this development are such that closed-loop configurations (generally 
referred to as nontopological) may be considered. If every body of the N-Body system is 
rigid, there may be no closed loops because such a system has an\indeterminate “load path.” 
To accommodate closed loops, the system must contain sufficient structural flexibility 
(compliance), and therefore modal displacement coordinates, so that the kinematic equations 
of interconnection constraint are algebraically consistent. 

The program development is such that none, several, or all bodies of the N-Body system 
may be flexible. The system may be “forced” by environmental factors such as gravity, 
gravity gradient, solar pressure, thermal gradient, and aerodynamic drag. 
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The computer program system described herein falls into several categories of capability: 

• Synthesis and time-domain solution of the nonlinear differential equations of motion 
of the complete spacecraft system idealized as a collection of interconnected flexible 
(or rigid) bodies 

• Linearization of the governing equations by numerical means 

• Time-domain solutions of the linearized equations that describe perturbation response 
of the complete spacecraft system about some predetermined (calculated or user- 
specified) nominal motion 

• General frequency-domain stability analysis corresponding to the linearized space- 
craft representation 

II. EQUATIONS OF STATE/TIME-DOMAIN SIMULATION 
A. Introductory Discussion 

The state equations* [that govern dynamic response of a system of interconnected flexible 
bodies, which may be actively or passively controlled and which may be “forced” by environ- 
mental factors such as solar pressure, gravity gradient, aerodynamic drag, etc,, are presented 


here in a concise summary form as: 

fill, - w;' ({G(i ♦ w/ N) «'■» 

Uli = PtijMi w-2) 

t^f ■ Z 1«, Ml (II-3) 

j 

l«l (iM) 

subject to the constraint equations 

E -M w-5) 

i 


♦For the interested reader, Reference Paper I describes the development of these equations in more detail. 
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In equations II-l through II-5 the index, j, ranges from 1 through the number of bodies of 
the system. Equations ll-l through II-4 represent n first order, nonlinear, ordinary differ- 
ential equations whereas equation II-S represents m additional conditions of kinematic con- 
straint. Thus, the dimension of the state space for a given system of controlled bodies is 
(n-m). That is, n-m state variables are required for defining the configuration at any instant 
of time, t 

State variables of the configuration space include absolute velocities, |u}j : modal displace- 
ments, |||j ; position coordinates (both angular and cartesian position), ; and additional 
variables,]^}- , that will subsequently be referred to as control variables. These variables are 
associated with the differential equations that define a given control law. However, they 
may reflect any other auxiliary differential equations that are necessary for defining the 
overall feedback system (for example, they may include thermal equilibrium states or other 
state variables necessary for a complete definition of the state-dependent environment). 

The right-hand sides of equations II- 1 through II-4 are functionally dependent on all the 
state variables and time, although the relationships may be termed only implicit at this point. 
Let it suffice that, in a way of introduction, a description of the nature of the governing 
equations II- 1 through II-5 be given here and more explicit development and discussion follow 
in subsequent sections. 

Equation II- 1 represents the dynamic equilibrium equations for the typical)^* body of the 
system. They are of the form shown whether the body is treated as rigid or flexible. They 
state, in effect, that a deformation dependent mass matrix, [m]., postmultiplied by a vector 
of relative accelerations, ful, , produces a vector of inertial forces that is balanced by all other 
state- and time-dependent forces, |G[., and interconnection constraint forces, [b]T |\|. The 
vector, jGJp includes inertial forces due to centrifugal and Coriolis acceleration, as well as 
elastic restoring forces, damping forces, control actuator forces, and so forth. The constraint 
forces, [b]T |x|, are necessary in order to satisfy the kinematic constraint equation (II-5); 
elements of the vector | Xj are actually Lagrange multipliers, evaluated and used in the solution 
process. 

Equation II-2 simply represents a selection transformation because the vector of modal velo- 
cities, ISli , is a subvector of |u}j. Equation II-3, used to develop represents a kine- 
matical transformation, transforming nonholonomic velocities to time derivatives of position 
coordinates. Finally, equation II-4 is an auxiliary differential equation that is user-defined 
and may be used to implement control dynamics and other feedback effects. 

The constraint equation (II-5) represents kinematic conditions of a form similar to those 
of equation II-3. In either case, there is a velocity transformation. Equation II-5 might be 
termed an active set of kinematic conditions, and those of equation II-3, a passive set. The 
active set is used to calculate m of the dependent elements of the |u[j vectors in terms of 
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the remaining independent elements and the prescribed velocities, | a|, some of which may 
be zero and others user-defined functions of time. Thus, the constraint equations are of a 
general form because nonholonomic, rheonomic conditions may be so represented. If the 
|u}j vectors satisfy the required conditions of equation II-5, the position rates, 1 13 }■, may be 
evaluated by the passive conditions of equation II-3, resulting in a kinematically consistent 
system. 

Note that m equations of constraint are represented by equation ITS. There are also m 
Lagrange multiphers in the vector, | X |. In studies of dynamic systems, the Lagrange multi- 
pliers and the dependent velocities and accelerations are often entirely eliminated from the 
governing equations. Such is not the case in this development in which Lagrange multipliers 
are involved in the equations for two reasons: (1) to monitor the multipliers as a function 
of system motion,, as they are interconnection forces and torques, and (2) for numerical 
implementation, it is convenient to calculate and use the { X [ vector in equation II-l . The 
Lagrange multipliers are calculated by differentiating equation ITS and combining that result 
with equation II-l giving: 

W ■ (E Ni ["’ll' wf) ' [{“} - E ( Sij Mi + Wi ["'ll' (II-6) 

Note the following functional dependencies: 

Wi-f (14, .{«},) 

m, ■ f(loMtli) (11-8) 

thus, 

Ih =f({4-l«l.M) (n-9) 


14, '((T1i) (ii-io) 
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( 11 - 11 ) 


Hi 


( 11 - 12 ) 




( 11 - 13 ) 


thus, 


[bi, 


( 11 - 14 ) 


{4 14 H- l«l’ 14 14^') 


( 11 - 15 ) 


and 


10} »f(lsf, idf, M- 1«(. 1?(. IH itj 


( 11 - 16 ) 


where, in the foregoing notation, the elements of the matrices/ vectors on the left are functions 
of the elements of the vectors on the right. The chronology of evaluations indicated must be 
followed in the solution process. 

The differential equations of motion for the system are therefore of the general form: 




( 11 - 17 ) 
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and the state vector and its time derivative are arranged as follows: 



where NB is the total number of bodies of the system, Nj3 is the total number of position 
coordinates necessary for orienting the system, and N6 is the total number of auxihafy 
(control) differential equations required. 

Now, if the |y| vector is known (numerically) from prescribed initial conditions or from 
numerical integration of | y[ , the primary task of the solution process is to numerically 
establish the |y[ vector. The |y[ vector is numerically (step by step) integrated to produce 
an incremented |y| vector, and thus a sequence of time-point solutions. 
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To summarize, a narrative description of the steps (numerical evaluations) necessary for pro- 
ducing |y^ given |y|, follows. 

The matrices, [BJ^ and [b]p are kinematic coefficients that depend on position and modal 
displacement variables and are evaluated as the first step. 

Now, if available numerical techniques (also computer software and hardware) were absolutely 
accurate, the |u|j vectors resulting from numerical integration of the vectors would 
satisfy the constraint equation II-5. Because this is not the case, the second step of the solu- 
tion process is to calculate the dependent elements of the |u|j vectors by using equation 
II-5. In fact, because of anticipated numerical inaccuracies, only the independent elements 
of the |u[j vectors are obtained by numerical integration. Only n-m “integrators” are 
involved in the solution process even though all of the elements of the |u[j vectors are 
numerically evaluated (by use of equation II-l); good numerical resolution is found in the 
independent |u}j elements because the Lagrange multipliers jX^were used. 

A kinematically consistent system results from satisfying equation II-5. The |U[j vectors 
may now be used with the selection and kinematic transformations as indicated by equations 
II-2 and II-3 to numerically produce all the modal velocities, and position coordinate 
rates, |^|, completing the third step of the process. 

Sufficient calculation has been completed to this point to evaluate the control variable rates 
according to equation II-4, producing |5|. During the process of calculating the |6|vector, 
all of the required control actuator torques (or forces) are calculated because sufficient 
numerical information is available. All of the constituents of the torques/force vectors, |G|j, 
are now available, and |G|., [m]j, and [b]j are therefore numerically evaluated, which 
completes the fourth step of the process. (Refer to the functional expressions of equations 
II- 1 1 through 11-14.) 

With reference to equation II-6, note that there is now sufficient numerical information to 
evaluate | X}, which is then used in equation II-l to calculate the |u}j, completing the fifth 
and final step of the process. 

Note that, in the above discussions, the solution process may be carried out through comple- 
tion, providing the state vector is numerically known. At any step of a simulation, the |y^ 
vector is known, of course, as the result of numerical integration. The initial state vector 
is another matter. It is difficult, if not impossible, for a user to prescribe |U[j vectors that 
are kinematically consistent with the conditions of equation II-5 ; also, the nonholonomic 
velocities of {u|j, when considered as a complete set, are of a somewhat abstract nature. 

The user is in a much better position to prescribe initial values of |j3[ and the initial 
velocities that are physically meaningful to him. Thus, to initiate the simulation (that is, to 
create an initial state vector from information the user is in a position to prescribe) the 
following preliminary steps must be taken. 
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• • 

The user must prescribe initial values of the ||3|, |l3|, and |S| vectors and the 

variable speed momentum- wheel spin velocities, 6 . Now, in that the prescribed position 
rates, |a|., are explicitly dependent on time and are always available, kinematic equations 
11-3 and II-5 may be used together to establish initial values of all |u[j. The question 
inevitably arises: Are the number of equations represented by II-3 and II-5 sufficient to 
solve for the elements of the juj^j? Consider the typical {Uj^j vector. Note that there are 
six reference-frame velocities in each {u (namely, cc , oj , co u, v, and w). Six relative 
velocities are also associated with each hinge. Now, if the system is a topological tree con- 
figuration, equations II-3 and II-5 comprise exactly the required number of equations to 
estabUsh the reference-frame velocities; that is, there are as many hinge points as there are 
bodies, and, even if every body were rigid, the system would be determinate. In this case, 
the initial sets of six reference-frame velocities are computed by equations II-3 and II-5 ; 
the prescribed initial ||| vectors and momentum- wheel spin velocities are simply placed 
in the appropriate ju} j vectors, and the initial state vector is thus defined. 

If the system is not a topological tree configuration, then there are more equations (II-3 and 
II-5) to be satisfied than there are reference frame velocities. (In other words, there are more 
hinges than bodies.) In this case, elements of the vectors must take on the responsibility 
of helping to satisfy the kinematic conditions. For each hinge in excess of the number of 
system bodies, there must be at least six deformation modes, represented by % coordinates, 
and they must be distributed throughout the system in such a way that the kinematic con- 
ditions of equation II-5 are independent. Qearly then, when there are more hinges than 
bodies (nontopological tree), one or more of the bodies must be flexible for the system to 
be determinate. When the configuration is nontopological, the user will specify initial values 
for all of the but he must acknowledge that they are not all independent and that the 
dependent ones (automatically determined by the program) are calculated and replace the 
values that he has specified. 

From these considerations, note that the initial state vector is created by the program from 
information that is user-supplied and that is physically meaningful to him. In the event of 
a nontopological tree configuration, the user’s only concern in regard to initial conditions is 
whether he has supplied an adequate description of system flexibility for the system’s 
kinematical equations to be determinate. 

B. Derivation of Equations of Dynamic Equilibrium 

The differential equations of motion and auxiliary equations that characterize a physical 
system may take any one of several equivalent forms. Equivalent form means that the same 
physical system can be characterized by more than one set of mathematical variables; in any 
case, the number of variables must be the same. For example, the motion equations for a 
rigid body could be derived by using Lagrange’s equations (resulting in six second-order 
equations), or the Newton-Euler equations could be used when translational motion is 
represented by three second-order equations, whereas rotational motion is represented by 
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six first-order equations (three moment-momentum equations and three attitude equations). 
In each case, there are twelve state variables. 

There are a variety of alternative methods of analytical dynamics that one may select from 
to develop the final (programmable) equation format A timely and valuable commentary 
accompanies the comprehensive comparative evaluation of these methods in a recent report 
by Likens (Reference 1). The basis for this development is effectively included in his 
discussion. 

The intent here is not to highlight any particular method of analytical dynamics as being 
superior to the others. Clearly, the methods are all equivalent if they are developed through 
completion without any compromising simplifications. The choice of method is made after 
considering the requirements associated with a particular problem or computer simulation 
program. This development begins with a Lagrangian approach ; then, algebraic manipula- 
tion is used to derive the format of equations Il-l through 11-5. 

Lagrange’s equations for the general situation appear as 


£ 

dt 



dT av V' 
— + — = Q. + > ^ 


Si \ 


for 0=1»2, • • n) 


(II-18a) 


n 

2 a., qj + a;, = 0 
j=i 

for (i= 1,2,* • in) 


(II-18b) 


In these equations, T and V are system kinetic and potential energies, respectively, and D is 
the Rayleigh dissipation function (accounting for internal damping). The generalized con- 
straint forces 


2 SiN 


augment the gener^ed forces, Qj (that arise because of the action of external factors), and 
are necessary for satisfying the additional conditions of constraint of equation II-18b. The 
form of equation 11-18 is complete and general, in that it includes unconservative forces 
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(explicitly time dependent), Qj, and dissipative forces, dD/dq^ and the auxiliary constraint 
equation (II- 18b) are in an all-encompassing form because holonomic conditions may be so 
represented. The coefficients, (a^, j = 1,2, •*, n; t), may depend explicitly on the time, (t); 
therefore, the constraint conditions as shown! account for both rheonomic and scleronomic 
situations. ' 

In the equations, n is the number of generalized coordinates involved in the representation, 
and m is the number of auxiliary conditions of constraint. Note that, although the q. are 
generalized coordinates (as they must be for the Lagrangian formulation), they are indepen- 
dent only in the isolated case when m = 0 or when there are no auxiliary constraint conditions. 
Some engineers share a misconception on this point: that, if the variables, qj, are not inde- 
pendent, they are not generalized coordinates. In view of the m constraint equations, this 
is simply a set of generalized coordinates that are not independent. 

In cases where all of the constraint equations are holonomic, it is theoretically possible to 
eliminate m of the qj in terms of the remaining n-m. However, if any of the constraint condi- 
tions are nonholonomic, a Lagrange multiplier, , must be used in conjunction with that 
equation. Of course, Lagrange multipliers may be used for either holonomic or nonholonomic 
constraints. 

In that the simulation program includes mathematical representation of active or passive con- 
trol for elements of the spacecraft system, some state equations involve control variables in 
addition to , equation 11-18. The manner in which the additional control equations enter into 
the composite system state equations is the same whether they are in the form given by 
equation II- 1 or that of equation 11-18. In either case, the control system state variables 
retain their identity, although the control forces/torques necessary to “close the loop” are 
transformed differently. In the case of Lagrange’s equations, the control torques contri- 
bute to the generalized forces, Qj, whereas, in the case of summary equation II-l , they con- 
tribute to elements of |G | and may be interpreted as ordinary forces or torques, acting at a 
structural hard point (or at a sensor point). Therefore, discussion of the control system 
will be postponed until later, and the “mainline” motion equations will be emphasized until 
the point is reached at which the control system coupUng can be clearly indicated. 

To “solve” Lagrange’s equations of motion, the exphcit form of the kinetic and potential 
energy functions, the dissipation function D, and the form of the transformation relating 
ordinary Cartesian position coordinates (positioning the typical system particle or element) 
to the generalized coordinates, qj must first be defined; the form of the transformation is 
necessary for expressing generalized forces, Q., in terms of external ordinary forces. After 

J 

the form of the energy functions and coordinate transformation is defined, the indicated 
differentiations (equation 11-18) are performed. A system of ordinary second-order differ- 
ential equations, which in many cases are nonlinear and which require solution using numeri- 
cal integration techniques, have been explicitly defined, but the motion equations have not 
yet been solved. 
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With numerical implementation and digital programming in mind, the form of the ordinary 
differential equations is recasted. First of all, they should result in canonical first-order 
form (the highest time derivatives appear uncoupled on the left-hand side). Also, compli- 
cated combinations of generalized velocities and displacements should be grouped so that 
such groups may be replaced with new variable names. These new variables have been 
called “quasi-coordinates” in the literature. This will simplify the required computer pro- 
gramming and minimize arithematic computation. Also, it helps considerably in organizing 
the numerical algorithms necessary for evaluating the left-hand side of the state equations. 
Thus, recasting the form of the governing equations is sufficiently justified. 

The recasting process is begun by defining the forms of kinetic and potential energy and the 
required transformation. First, note that bodies of the system of flexible bodies are tenta- 
tively treated as if they are completely independent of each other. The influence of any 
body on another is accounted for by the additional constraint conditions and the Lagrange 
multiphers. Thus, if kinetic and potential energies for the typical body are expressed and 
Lagrange’s equations are applied to it, the ordinary differential equations pertaining to it 
are simply a subset of equation, II- 18, and the total system through the representative form 
of the typical body will have been accounted for. 

The generalized coordinates chosen to represent the configuration of the typical body include 
three Euler angles to indicate attitude of the body fixed-axis system relative to an inertial 
frame, three projections (components) of the position vector from the origin of the inertial 
frame to the origin of the body fixed-reference system onto the inertial axes, and N elastic 
displacement coordinates. Note that the origin of the body fixed-axis system need not 
necessarily coincide with the body center of mass. Also, the elastic displacement coordinates 
may be measurements of displacement at a discrete set of points on the body, or they may 
be coordinates associated with normal vibration modes. In either case, they represent dis- 
placements measured in the body axis system. For the flexible body, its generalized coori- 
dinates are tabulated as: 



Attitude 
Euler Angles 


Body’s Reference Point 
Position Coordinates 


^1 

^2 


Elastic Displacement 
Coordinates 
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A transformation now exists that relates a set of nonholonomic velocities to the generalized 
velocities that are extensively used in recasting the equations. The transformation appears 
as; 
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(11-19) 


where, in equation 11-19, the vector of nonholonomic velocities, |u }, contains the three 
projections, (co^^, cu^), of the angular velocity vector, co, onto the body fixed axes 

(cj is the angular velocity of the body reference frame), the three projections of the refer- 
ence point translational velocity, (u, v, w), onto the body fixed axes, and the displacement 
rates, | ^ The elements of the transformation, 7 jj (i, j = 1,2, 3), are direction cosines; 
the submatrix, [ 7 ] , is an orthonormal rotation transformation relating the attitude of the 
body fixed-axis system to the inertial frame. The submatrix, [tt] , is also a rotation transfor- 
mation; however, it is not orthonormal because it relates vector components based on an 
orthogonal basis to those of a skew (nonorthogonal) basis (namely, the axes about which 
Euler rotations are measured). 

In short. 


w ii) 


( 11 - 20 ) 
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Clearly the elements of [jS] are functions of the three Euler angles. Twelve sets of Euler 
angles are possible. Any one set is valid for use in subsequent development; the resulting 
equation form is independent of selection from the twelve sets of angles. 

Elements of the transformation, [/3], may be explicitly defined in terms of three of the 
generalized coordinates (the Euler angles). 

The kinetic-energy expression for the body is most easily expressed (initially) in terms 
of the nonholonomic velocities, |u[, after which [/?] is used to replace |u[ with [|3] |q^. 
The kinetic energy is then expressed completely in terms of generalized displacements and 
velocities— the form necessary for applying equation 11-18. 

Kinetic energy for the typical body is 


T = 



V vadV 


( 11 - 21 ) 


where v is the velocity field, and a is mass density, and where integration is carried out over 
the volume, V, of the body. 

The inertial position of any point, p, of the body (figure 3) is: 


r 


Xr 


+ Po + n 


( 11 - 22 ) 


where Xj^ is the inertial position of the body’s reference point (R is the origin of the body 
axis system), positions the point p' (which coincides with p in the undeformed configura- 
tion) from point R, and r] (x, y, z, t) is a measure of elastic displacement. 

The vectors Pg and rj are referenced to the body axis system, thus 




(11-23) 
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Figure 3. The flexible body. 


and 


N. 

n (x, y, z, t) = I T T k J / 

k=l ' 


V y> 2) 

,^zk y> 2)J 


fk (0 


( 11 - 24 ) 


the elastic displacement rj is represented as the superposition of a finite number of single 
valued space functions, 

The velocity field, v, is obtained as 


= '^R + (^0 + ’j ) + 13 ^ 

k=l 


( 11 - 25 ) 


with 


- _ ^ 
dt 


18 



The velocity of the reference point, R, may be expressed in terms of components referenced 
to either the inertial frame or the body frame; that is, 


also 



(11-26) 



V 

w 


The unit vectors, I i , j, k}, {l, J, K }, are related through the rotation transformation, [ 7 ] , 
and it follows that 


*u" 


"X" 

V 

= [ 7 ] 

• 

Y 



• 

w 


_Z_ 


(11-27) 


To be concise, the repeated index summation convention is introduced at this point. With 
this convention, when any two factors of a term have the same index, summation over the 
range of that index is implied, and the 2 sign is deleted. For example, the third term on the 
right of equation 11-25 is 


and represents 

thk 

k=l 
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Now, if equation 11-25 is substituted into equation 11-21, the kinetic energy is 




T = I ^ + [<o X (p^ + T?)] • [OJ X (p^ + T?)] 

V 

+ 0k • ^ ik 

+ 2 V ^ • [o X (pg + ij)] + 2v^ • 0^ 

+ 2 [w X (pg + ij")] * 0k Ik 1 


or, integrating term by term over V, 


T = V4mj^uvwJ |uvw| 
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(11-28) 


(11-29) 
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+ 



(II-29) continued 




m 



(11-30) 
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(11-31) 
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and 



(11-32) 
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Also used was 


f' 


= / ^Xk 


(11-33) 


and 


Bjk =J ^xk + ^k + ^zkj 


(11-34) 


f (x + S. ) adV 


(11-35) 


= S,o + a^. 


and also, 


•■xk 


|Sjzk ' ®zjyk j 


adV 


^yzk ^zyk ^ I ^yjzk * ^zjyk 


f ^ 


xy 


J +(b . + b .|£+c..t£, 

xyo y xyj yxj j xjyk 


(11-36) 


(11-37) 


All other quantities involved in equation 11-29 are obtained by cyclic permutation of the 
indexes, x, y, and z. Finally, because the kinetic energy is of quadratic form in the elements 
of |u }, it may be expressed as a triple matrix product 


T = % [uj [m] |U} 


(11-38) 
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or in short, 
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Using equations 11-40, 11-19, and 11-38 gives 

T= \q\wf M m Jq( 


(11-39) 


(11-40) 


(11-41) 


Clearly, the elements of [m] depend only on the ; the elements of [J3] depend on the Euler 
angles, and kinetic energy is therefore a function of generalized velocities and the generalized 
coordinates themselves. Thus, the functional notation, 

T = T (qi.q2.-*q„;q,,q2,-*q„j 
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is applicable; terms such as 9T/3q^ will come about and play an important role in the 
simulation. 

To continue, it is necessary to express the potential energy, V, and dissipation function, D. 
Assume that the elastic strain energy can be written as a positive definite quadratic form in 
the elastic displacement coordinates, or 

V = [fj [k] j (11-42) 

the symmetric matrix, [k] , is developed by standard finite-element techniques such as those 
embodied in NASTRAN. If } is a set of normal modal coordinates, then [k] is diagonal 
with the diagonal element appearing as 


k„ = (11-43) 

where co. is the natural frequency. Of course, normalization of the eigenvectors (mode 
shapes) is assumed so that the generalized mass for the vibration mode is unity. 

Now, because 


l«f= [0j0|l„l (q( 

■ iStHi} 


(11-44) 


it follows that 


y = 'A LqJ[Sil^ [k] [Sj] {q( (11-45) 

Similarly, D is written as 

D = [C] [Sj] (11-46) 


where matrix, [C] , developed by standard finite-element techniques, is equivalent viscous 
damping for the structure. 
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Refer back to Lagrange’s equations (II-18a and II-18b), and reexpress them in matrix format 
as follows: 


H M|q()=-[S^l^ ([k] lSj]{qf+ [C] [S^])q}) 

+ |q|+ y]^ [m] M|q}l 

(IM7) 

+ */^{[qJ I^]’’ H [/3,j]|q}} 

+ ^ {[qj [ni j] Mlqlj+^faj'^jx} 

and 

[a]{q}=-|a,| (11-48) 

What is meant by [^ , ] and [m , ] is the partial derivative of each element of [|3] and [m] 
with respect to the jf" generalized coordinate. 

Now, define the ordinary momenta (see Reference Paper II); 


|p}= H [)3] |q} 

= [m] |u} 


(11-49) 


Also, because 

M= m |q( 


it follows that 


tq(= W Wl 


(11-50) 
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Using equations 11-49, 11-50, 11-47, and 11-48, 


{pS= [Sj]jq|+ [C] [Sjljqf) 

+ H M |LuJ KiJ M}^ 

[a] [^] 'M = {-a,l (n-52) 


On studying equations II-5 1 and 11-52, several observations can be made: First of all, recall 
the form of 1/3] and [S^ ] (equations 11-19 and 11-44). It is clear from these forms that 


[/3]-»t [Sj]T ^ 


(11-53) 


and that 


[S^]|q}=|l} (11-54) 

and 

[Sj] |qf=li} (11-55) 


Also, because the elements of [m] depend only on the first six elements of 

||_uj (m jl )ul| 


are null; therefore 


HI " {LuJ ["■.)! ju(} “ (LuJ [nij )u|} 


(11-56) 
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Further, note that matrix transforms the generalized forces, |QJ, to forces “acting in 

the quasi-coordinates,” or call 

|o„l - wr"iQ(, (11-57) 

thus, jGgjjj- contains ordinary forces and moments attributable to external sources and cor- 
responds to time derivatives of the ordinary momenta. 

Because the transformation, [/J] , depends only on the Euler angles, it follows that only the 
first six elements of the column 

("tj- 

are nonzero, and, after considerable algebraic manipulation, this column may be reexpressed 
as 

[H] 
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With these observations and definitions, equations 11-51 and 11-52 may be reexpressed as 


= KJ - [“-] m - [e] in - M H 

+ ^ |LUJ [m.j] {uf} + [b]^|X[ 


and 


w jut = in 


where 


[b] = [a] [p]-‘ 


was used, and 


{“} =-w 


(11-59) 


(11-60) 


(11-61) 


(11-62) 


Note that constraint equation 11-60 is now expressed in terms of the nonholonomic velocities, 
{u[ ; the coefficients, [b], are obtained directly from relatively simple, vectorial expressions 
of kinematic constraint. The same [b] coefficients are transposed and used to multiply 
{x[, producing constraint forces/torques corresponding to the ordinary momenta. 

If the I G f vector is now defined as 


1<=! = 1g„( - li} - [j] li( + (5) H )u[ 
* « IluJIiHj) 1u(| - |ii)'lul 


(11-63) 


it follows that dynamic equilibrium equations for the typical r^* body may be written as 


i 








(11-64) 
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to be used in conjunction with system kinematic constraint equations 


^[b|, {uj,. {«( (11-65) 


which is the same form as that given by equations II- 1 and II-5. 

The last three terms of |G| given in equation 11-63 are inertial forces that involve velocities 
and displacements of the body. The matrix [m] is an instantaneous inertia matrix, depend- 
ing on instantaneous values of the deformation coordinates, | ^ The centrifugal and 
Coriolis effects are completely accounted for within the framework of the assumed velocity 
field (given by equation 11-25). These effects would not be accounted for if “tangential” 
velocity due to elastic displacement were neglected (i.e., if it were assumed that 
I CO X I « I 6j X Pg I. In this case, the inertia would be constant and independent of 

An accurate definition of the dynamic equilibrium equations clearly hinges on a complete 
and accurate definition of the constituents of the vector, which includes the inertia 
matrix, [m]^. Also, the kinematic coefficients, [b]^, must be developed in an exact fashion. 
Kinematics and a more explicit development of |g} are given in subsequent sections. 

C. Kinematics and System Topology 

From a Lagrangian formulation, all of the generalized forces, not derivable from a potential 
function, ordinarily appear as |Q | on the right side of Lagrange’s equations of motion. 
Internal damping forces have been accounted for with the use of Rayleigh’s dissipation func- 
tion, D, and, for generalized constraint forces, by Lagrange’s multipliers. 

Thus, the remaining generalized forces to deal with include those that are attributable to 
external factors such as aerodynamic drag, solar pressure, and other commonly encountered 
environmental loadings. 

Control forces (servodrive torques, reaction jets, etc.) are also treated as if they are external. 
They are not explicitly external, of course, because they depend on time through position 
and rate errors that are functions of elements of the state vector and on control system 
state variables that arise from a given control law. 

Assume that there is a finite number of points on the typical body at which a force vector 
(or torque) is known to act. Each of these force/torque vectors contributes to the generalized 
forces, |q|. The generalized forces are calculated by expressing the virtual work of the ex- 
ternal ordinary forces in terms of virtual displacements of the points of force application. 

The transformation that relates ordinary coordinates to generalized coordinates is then used 
to define the explicit form of the generalized forces. For example, suppose that a force fp, 
and torque, Tp , act at point p of the typical body. Their virtual work is 

6W = f • St + T *60 (11-66) 

p p p p 
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Note that the virtual rotation, 8^^, was treated as a vector quantity. This is valid, even though 
a general rotation is not a vector quantity, because the virtual rotation is infinitesimal and 
is therefore a vector. Further, because virtual displacements are infinitesimal, 6 and 5^p 
may be expressed in terms of virtual displacements of the quasi-coordinates; that is. 


6r_ = 


b j kj ^ 
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.*^3. 
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-(z + 7} ) 0 

p *zp' 
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(11-67) 


and 


89. 


-Lr r iy( 


~89' 

X 

+ 


50, 



60 

Z 


%(Xp.yp.Zp) 




) 


(11-68) 


where (5r^, Spj, Srj) are components of virtual displacement of the body’s reference point, 
R; (S0.,, 50„, 60,) are components of virtual rotation of the body axis system; and {a., 
a., a.) are components of the space function, a., representing elastic rotation at point, 

yj zj j 

p (modal slopes, for example). 

Now, assume that the force and torque vectors, (fp and Tp), are referenced to the body 
axis system; they may therefore be written as 



(11-69) 
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and 



, Note that virtual displacements of the quasi-coordinates are related to virtual generalized 

displacements by the same transformation that relates nonholonomic velocities to generalized 

velocities (equation 11-19). It follows that the virtual work attributable to f and T maybe 

PP 

wntten as 



The virtual work is also expressed as 

6W = [Sqj |q[ 
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and, because Sq^ is arbitrary and independent (treated as independent in the face of Lagrange 
multipliers and constraint equations), it follows that 


|q} = M'" [bpl’’ 



(11-72) 


Equations 11-71 or 11-72 have a noteworthy geometrical interpretation. Note that the first 
three lines of 




T 



are components of the resultant torque vector, Tp + (p^ + rj) X fp, acting at the body’s 
reference point, R. The second three lines are components of the resultant force vector, fp , 
whereas the line, 0 > 6), corresponds to the standard procedure (of structural dynam- 
icists) for calculating Q^, or, as it is usually expressed, generalized forces acting in deforma- 
tion modes are 

|q} =[«i’'1t} 


Also, recalling the form of [/3] (equation 11-19), note that [tt]^ resolves the resultant torque 
vector (about orthogonal body axes) to components about skew axes about which Euler 
rotations are measured, whereas [7] ^ resolves the resultant force vector (about orthogonal 
body axes) to components along the inertial axes. Further, note that [bp ] is a matrix of 
coefficients that relates the velocity of any point, p, to the vector, |u[. This provides ad- 
ditional insight as to why the same coefficients that are used in kinematic constraint equa- 
tion 11-60 are used (in transposed form) to multiply -producing resultant constraint 
forces. 

Thus, the remarkable duality of purpose associated with [b]-type coefficients has been 
emphasized. They are initially expressed by writing simple kinematic velocity relationships. 
The coefficients, [b] are then used to transform discrete ordinary forces and torques to 
equivalent forces and torques acting throu^ the body’s reference point, R. The matrix, 

[|3] , which is also a velocity transformation, is transposed to produce the transformation 
to generalized forces (if they are desired). 
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For ordinary momenta equations, the desire is simply to express | | , which (following 

equation 11-57) is given by 



(11-73) 


This p given by equation 11-73 reflects only the contribution of the force/ torque acting 
at a sin^e point, p. The total must be obtained by summing all the points of the body 

at which forces and torques act, or 


NP 

1=1 


V 


(11-74) 


Kinematic coefficients, [bp ] , such as those of the previous example, will be required through- 
out in the formulation of the state equations. They are used to synthesize the constraint 
equations and to produce |Gj^, and they are involved in the velocity transformation of equation 
II-3. It is therefore advantageous to think of a “bank” or collection of all the required kine- 
matic coefficients to be put together in a semiautomatic fashion by using input specifications 
to the digital program. 

7. Sensor Point Kinematics— Force/Torque Transformations 

Consider the typical structural hard point, s (figure 4). Assume that a right-handed triad is 
fixed to point s and that the elements of the triad are unit vectors labeled 1, m, and n. Now, 
body, n (which has point s on it), also has a right-handed triad fixed to point n. Suppose 
that, even when body n is in an undeformed state, the s-triad is misaligned with respect to 
the n-triad. When the body deforms, there may be further angular misalignment between 
the two triads. Thus, the relationship linking the two sets of unit vectors is 





(11-75) 
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Figure 4. Two typical contiguous bodies of the system. 

where [jR^i] and are orthonormal rotation transformations, the first relating the 

“naturally” misaligned triads via constant Euler rotations and the second accounting for 
additional rotation due to the body’s deformation point s. 

The structural deformation at point s is assumed to be sufficiently small that the Euler 
rotations associated with may be evaluated through the use of 



Ki id 


(11-76) 


where [aj is a (3 X N) matrix of modal rotation amplitudes at point s. (Each of the N 
columns corresponds to a deformation mode.) Concisely' denote the traids associated with 
points n and s by { e^^ | and { [, respectively. The relationship Unking the two sets of unit 

vectors may then be expressed as 


f'.t - [.R„] I'd ai-77) 

In subsequent kinematic development, there is a requirement for expressing the absolute 
velocity of a typical s-point and the angular velocity of the- typical s-triad in terms of velocity 
states of a given body. Picture a six-long vector (column) of velocity components (three 
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rotational and three translational) that are projections of and onto the s-traid axes. It 
is related to the {U|jj vector for the body by the transformation 
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(11-78) 



where [hj and [a^] represent matrices of displacement and rotation amplitudes, respectively, 
and [S(">] is an antisymmetric matrix accounting for a vector cross product, or 
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(11-79) 


The superscripts in equations 11-78 and 11-79 are used to indicate the frame to which the 
velocity components are referenced. 

Kinematic coefficients such as those of equation 11-78 are generated for each so-called sensor 
point of the system of bodies. They are used by the simulation program to produce contri- 
butions to |Gg^ I from given force/torque components in the manner indicated by equation 
11-74. 
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2. Hinge-Point Kinematics 

Kinematics associated with hinges follows a line of development somewhat similar to that of 
sensor points. Consider the points, p and q (figure 4), to be two structural hard points 
associated with a given hinge. All necessary kinematic information pertinent to the hinge 
is obtained through expressing the velocity of point q relative to point p and in expressing 
the relative angular velocity between the q- and p-frames. It is convenient that the angular 
velocity components are projections onto skew axes (Euler angle rates) and that translational 
I velocity components are projections onto the axes of the p-triad. The six relative velocity 
i components may be assembled into a column matrix as 



(11-80) 


where represents the three relative Euler angle rates and represents the three 
relative translational velocity components all of which pertain to the k^* hinge. The column 
of relative velocities may now be expressed as 


(b,l, lu|. 


(11-81) 


with 


[bp] = 


W" [.M [01 I - [P^m] [%] 


- [p^nG [SKl. '-[pRj -[pRnG [bp] 


(11-82) 


and 


[b,] = 


[.^n] 


[ 0 ] 


[pM m [sS]|[p^] m 




[pM [.^n] [bj 


(11-83) 


In equations 11-82 and 11-83, the rotation transformations, [pR^ ] and [qR„l , are developed 
to include the effects of structural deformation as indicated in equation 11-75; the rotation 
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transformations, and [pRq 1 . are developed in standard fashion using the three Euler 
rotations, 

For further discussion, consider the system of bodies shown in figure 5. Topology of the 
I system is simply indicated by an integer array, called “ITOPOL,” as follows: 


[ITOPOL] = 


12345678 — Hinge number 

— Body (n) relative to 
— Body (m) 


1 2 4 3 5 6 7 7 
0 3 2 6 3 1 5 2 


The [ITOPOL] array, which is the actual input to the simulation program, is used to define 
system topology as indicated. Now, with reference to the example shown in figure 5 and 
the corresponding (ITOPOL) array, the form of the velocity transformation may be written 



|U)4 

M, 


(11-84) 

We 

W, 

Ws 
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NOTE: 


Hinge labels are circled; 
body labels are not circled. 



where [bpj j] and [bqj j] are matrices as defined in equations 11-82 and 11-83 (with i = 

Hinge number and j = Body number). The velocity tr^sformations of equation 11-84 repre- 
sent the “bank” of all hinge kinematic coefficients previously mentioned and produces every 
possible velocity component pertinent to hinges. Referring to basic system equations II-3 
and II-5, note that selected lines, or equations, from the bank (equation 11-84) are taken to 
represent constraint equations or position coordinate rate equations. The [B]j and [b]j ! 
coefficients of equations 11-3 and II-5 are simply subpartitions extracted from equation 11-84. 

To implement calculation of Lagrange’s multipliers (equation II-6), it is necessary to develop 
time derivatives of [b] , coefficients. In a manner similar to the foregoing, in which all [b] . 

J • J 

coefficients are extracted from the complete collections, the [b], matrices come from a 

• • J 

collection of matrices whose members are [bq. .] and [bp. .], which are developed in 
Appendix C. 

D. Development of the {g[. Force Vector 

The equations of dynamic equilibrium for the body of the system are given in an earlier 
section as equation II- 1. As noted there, the right-hand side includes a so-called {cfj vector, 
which accounts for all state-dependent forces except for those of interconnection constaint. 
In equation 11-63, the vector is presented in a somewhat more-developed form. 
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The purpose of this section is to provide more explicit development of the elements that 
contribute to {Glj. All contributions in the following expression may be accounted for 
(omitting the j subscript, understanding that the typical, or j'* body, is being dealt with): 





+ [S2] [m] 



+ 




- M + 1g„J + jG„( 


(11-85) 


Although the first term,|G^^| , has been discussed in the previous section (equation 11-74), 
note here that the ordinary force/torque components that produce |Ggj^[ may be considered 
as a miscellaneous force vector. Its presence provides the program user latitude for including 
a variety of additional effects. Clearly, it is the implement through which control forces/ 
torques are “fed back” to the dynamic system. 

The second and third terms of equation 11-85 have been previously introduced. There is no 
implicit restriction on the stiffness and damping matrices, [k] and [C] ,^nor is there a restric- 
tion on definition of the 1^1 coordinates; in the majority of cases, they will likely be coordi- 
nates associated with orthonormal vibration modes. However, they may be physical (ordinary- 
discrete) displacement coordinates as well. In the latter case, the [k] and [C] matrices are 
usually coupled. 

The last two terms of equation 11-85 are included to account for momentum-wheel coupling 
and gravity effects, respectively. The treatment given to built-in momentum wheels is such 
that, in addition to producing a contribution to |G^, there is also a required extension to 
the form of the [m]j matrices because momentum wheels are inertially coupled. Thus, there 
is sufficient requirement for a dedicated development concerning momentum wheels. The 
next two sections deal exclusively with momentum-wheel and gravity effects, respectively. 

The remaining terms that contribute to |G[ are basic inertial effects and involve the matrices, 
[m] , [ni ^ ] , and [in] . With reference to equation 11-39, the form, [m] , is given, correspond- 
ing to the case with available single-valued space functions 0^. Ordinarily, access to such a 
description of the structure’s deformation modes is not possible because of the structural 
complexity of typical spacecraft. The analyst should always be able to obtain, as data, 
matrices of mbdal amplitude ratios (mode shapes) and the corresponding structural mass 
matrix (generated by finite element techniques). To accommodate data based on the more 
practical definition of structural characteristics, it is necessary to recast the inertia matrices, 
[m] , in a similar but more general format The generality of the development of paragraph 
II.B is not compromised by extending the form of the inertia matrix. The extended, or more 
general, inertia matrix is developed in Appendix A, but here, for purposes of developing 
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inertial contributions to the |g| vector, the resulting form is accepted and the kinetic energy 
expression is presented as 


T = 14 [u] ([mj + {v\ (11-86) 


with the repeated index summation convention implied, and with [m^ ] of the form 


KJ 


I -S I d 
1 m I a 

jT ' T 1 
[d ,a , e_| 


(11-87) 


that is, it is identical to the [m] given by equation 11-39 except that it is constant, independent 
of deformation. The constant inertia matrix, [m^ ] , as given by equation 11-87, is always of 
the form shown regardless of the choice of “modal” columns. The form of the matrices, 

[m^ ] and [m 2 ] , is such as to accommodate the general situation; that is, their definition in- 
cludes inertial integrals as defined for a continuous system (equations IT30 through 11-37) 
or as defined by structural mass matrices that are called “lumped” or “consistent.” 

The inertia matrix associated with is 


2 b, -b, -b. 

<k tt a 

1 2 3 



“4 “s “6 

IfzAil 


“7 “8 “9 



0 0 0 

000 


0 0 

000 


0 

000 



000 

(symmetric) 


0 0 



0 





( 11 - 88 ) 
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and the one associated with is 


Kljk = 


C -C -C 



11 12 13 



C -C 

'“22 '^23 

0 

0 

c 
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0 

0 

(symmetric) 


0 


(11-89) 


Now, for N-deformation modes associated with a given body, it is understood that the range 
of the indices, j and k, is N; the coefficients, *• therefore stored 

as 9 (N X N) arrays of inertial integrals, whereas (bj)., C^)., •• (bgX and («j )., ( 02 )., *• (a^). 
are stored as a (6 X N) array and a (9 X N) array, respectively. Thus^, from a programming 
standpoint, note that 9N^ + 1 5N storage locations are required for accommodating the inertial 
integrals necessary to account for the deformation-dependent mass matrix. Of course, if a 
particular body is rigid (N = 0), then only the first (6 X 6) diagonal partition of [m^ ] is used. 

When the body is flexible (N >'0), the inertia matrix is calculated from deformation states 
(^j) and inertia integrals in the manner indicated by equation 11-86; the redundant operations 
due to symmetry and null operations are avoided in the digital code. 

Having an instantaneous numerical evaluation of the inertia matrix, the term, [f2] [m] |u[, 
is calculated and added to |g[, consistent with the expression of equation 11-58. 

It is now possible to express explicitly, the combination of the remaining two inertial force 
vectors in terms of the inertial integrals given in equations 11-88 and 11-89. For further 
development, the combination may be defined as 

IgJ = I LUj |U(| 


[■ill {u( 


(11-90) 


41 



Thus, the first element of |G^.| , corresponding to is 

(G.), = I -2<^, (b,), * (b,), + (b,), 

-u(ttj)j - - w(a3>j 

( 11 - 91 ) 

■*'Wy l(^i2)fij ■*■ (^12^^ ^ ‘^z ^ ^8 | 

The second element, corresponding to co ^ , is 

= I (^)j -^‘^y (ba)j + (l>6)j 

-u(%)j - v(aj)j - w(ag)j 

( 11 - 92 ) 

■(^zxV^k ^x ^8 

-2Wy (^22)5] ^8 ■'■ t(^23^8j '*’ (^23^^ ^8 I 

The third element, corresponding to <0^ , is 

(Gc)3 = I <^x (^s)j + S (be)j - 2Wz (^3>j 

-u(oi,)j -v(otg)j -w(a^)j 

( 11 - 93 ) 

xyV + Wjj K^13^8j ^ ^8 

+ [(^23^8} '*’ (^23^J8l ^8 ^^33^8j ^8 | 


42 



The fourth element, corresponding to u, is 


(G,)^ = + Wy (a^)j + (a,)j | 

The fifth element, corresponding to v, is 

(Gc)j = - (o(2)j + cOy (aj)j + (ttg)j I ki 

The sixth element, corresponding to w, is 

(Gc) 6 = («3)j + cOy + co^ (a^)j j 

Finally, for the element k + 6, corresponding to an inertial force acting in the 

+ "J KC„)y«, ♦ (VJ 

K^33^kj ^ 

~ I KG23)fcj + (G23)jkl ^ ('’e^kj 

+ [(a,)^ u + V + ia^\ w] 


( 11 - 94 ) 

( 11 - 95 ) 

( 11 - 96 ) 

coordinate, 

( 11 - 97 ) 
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(11-97) continued 


+ cOy K%\ u + V + (a^\ w] 

+ [(a,)^ u + (a^\ V + (a,)^ w] 

+ K^yz^kj " ^ ‘^y K^zx^j ■ 

+ “z KC,,)y - (C«,),J } ( 


On examining the composition of the inertial force (G(.)^^.g, note that the first six bracketed 
terms represent centrifugal forces (distance X omega-squared) acting in the deformation co- 
ordinates, whereas the last bracketed terms of equation 11-97 represent Coriolis forces 
(velocity X omega). 

E. Momentum-Wheel Coupling 

The spacecraft system undergoing analysis may have several “built-in” momentum wheels. 

A momentum wheel is usually defined as a cylindrical or disk-shaped mass that spins about 
an axis that is fixed to a structural hard point of a given body. The wheel can be either 
spun up or despun by an electric motor whose rotor is part of the rotating mass. The shaft 
torque that acts to accelerate the wheel also acts on the body in a negative sense, providing 
active attitude control. The shaft torque is generally governed by a control law th^t “senses” 
attitude and rate errors of the body. In this development, a momentum wheel is assumed 
to be inertially symmetric about its spin axis. 

To develop the inertial coupling effects of the typical momentum wheel, consider three 
unit-vector bases: 


Le„J = [ i, j, kj 

(11-98) 

LcjJ = [e.m.nj 

(11-99) 

Le^J = n'J 

(II- 100) 
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The first triad is the body-reference triad for body n, the second is a sensor-point triad (fixed 
to point s), and the third triad is fixed in the momentum wheel. One of the three unit vec- 
tors of (_e^J is coincident with one of the unit vectors of either C, m, or n may be 

the spin axis depending on the preference of the analyst). In figure 6, n = n' is shown as 
the common, or spin, axis. 

The absolute angular velocity of the L®wJ fr^me can be expressed as 

-- W (i“4 ^ l^l ») 

where is an elementary three-long position vector (null except for unity in the first, 
second, or third locations corresponding to fi , m, or n being the spin axis), and 0 is the 
relative angular speed of the with respect to the L^gJ frame. 

With the inertial characteristics assumed (axisymmetry) for the wheel and with the velocity 
expression of equation II-lOl, the total angular momentum vector for the wheel may be 
written as 


!> = Uw) l"wl 

( 11 - 102 ) 

- I'J I'wl ( t“ f + {Pwl «) 

with [j^] diagonal with all diagonal values equal to Jj except for the position corresponding 
to the spin axis, J^. J-j. is the mass moment of inertia about any axis perpendicular to the 
spin axis, and is the spin inertia for the wheel. 



Figure 6. Typical body/momentum-wheel relationship. 
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The torque acting on the wheel (resolved to the [ e^ ] frame) is 


T ' LnJ N - ^ 


( 11 - 103 ) 


= Le.J ( [J„I I^J * jP„f !,«' - [fi.i (J„I Kl - [£2.1 {Pj J.») 


where an SK* operator is defined so that 


[J2J = SK* {coj 


or 




— 


— — 

0 

^S3 

—Cl) 

S2 



-6J 

^S3 

0 

sl 

= SK* 


CJ 

^s2 


0 









( 11 - 104 ) 


The torque acting on body n at point s due to the wheel is -T, and it drives the body’s quasi- 
coordinate as 


• - Tb.]" (dwl fb.l (ui„ 

- (n.i [J„l {<p. I - HI 


* HI [b.l K + 

Kl '/) 


Ip i j 0 

l w> s 


( 11 - 105 ) 


with 


A] 

■W [’ir-] 

(11-106) 

and also, as can easily be shown, 



(b.l |U} . = 

(sk*([.R„] KI lif„)l [b.l |u}„ 

(11-107) 
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Now, the shaft torque is simply the projection of T onto the spin axis, or 


\ ‘ [>■»] {Tf 

■ [PwJ ItJ iu|„ * J. [P„J lb.) )u}„ * 


(11-108) 


Equations 11-105 and 11-108 permit the coupled equations for body n and several momentum 
wheels to be expressed as 



{G}„ + [b]^ |X} 


X Uwi Kl - S fb,r Uwi h 

w w 

0 

+ 

0 

0 


0 


(11-109) 


+ 


-X iyT(SK*{pj)rb,] 


- Js. [Pw.J Ib.l 



The inertially coupled body/momentum-wheel equations (for two wheels) are shown as 
equation 11-109 simply for the purpose of indicating the form. Note that, within the 
equations, there effectively resides the original form of the dynamic-equilibrium equations 
for body n; namely. 


M„ {u(„ = (g|„ + [h]l {x} (ii-iio) 


which govern if no momentum wheels are associated with body n. In equation II-l 10, the 
caret (A) has been placed over G to represent the right-hand side force vector that excludes 
momentum-wheel effects. 
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Now, on further study of the form of equation 11-109, note that, if the “locked” momentum- 
wheel effects are already included in the definition of [ni] (which is the standard practice 
when inertially coupling systems together), the (1, 1) partition of the coefficients on the left 
of equation 11-109 becomes simply [m]„. Also, the second column on the right of equation 
11-109 is absorbed in having already been accounted for in the development of dynamic- 

equihbrium equations. 

It therefore follows that, in order to implement momentum-wheel coupling with one of the 
flexible bodies, it is only necessary to extend the vector to contain momentum- wheel 

spin values, (0), to extend the inertia (except for the [ 1, 1 ] partition) as indicated in equation 
11-109 and to add to the right-hand side force vector 



W 

T.. -J.. luj 
1'.. - ’.2 [PwaJ luj 


(IM 1 1) 


The values for shaft torque that appear in are established by a given control law 

if the wheels are to be considered variable speed. If a given momentum wheel is of constant 
speed (used only for “gyroscopic damping”), the torque equation for it is deleted from 
the form of equation 11-109; however, its effects are still included in the upper partition of 
the vector, (the gyroscopic torque due to constant 0 ). 

Clearly, the equations of dynamic equilibrium for a body, after having been augmented to 
include momentum-wheel coupling, are still of the general form 

|u|. = [m]:‘ + [h]J (11-112) 


F. Gravity-Gradient Effects 

Attitude dynamics of orbiting spacecraft can be significantly influenced by the gravitational 
force that is distributed according to the system’s position and deformation state. The 
gravitational force per unit mass varies (in a central force field) simply because different 
mass particles are at different distances from the Earth’s center of mass. Figure 7 describes 
the geometry associated with a typical elastic body. 
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Figure 7. Geometry for gravity effects on a typical body. 
For a central force field, the gravitational force per unit mass is given as 



which, to a first-order approximation, is 


( 11 - 113 ) 



wherei 

GM = the Earth’s gravitational constant 
mj = the typical mass particle 
g^ = local gravitational acceleration 
e"„ = a unit vector directed along R„ 

c = the origin of the body reference system 
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The virtual work due to gravitational force can be written as 


5W 


g 




(11-115) 


with nij replaced by differential mass adV. 

The virtual displacement field is expressed in terms of virtual displacements of the quasi- 
coordinates as 


5r = 5r^ + X (p^ + t?) + 6t? 


(11-116) 


In combining equation II-l 15 with equation II-l 16, the torque about point c, due to gravityr 
gradient effects, is 



(11-117) 


where 

S = the first mass moment about point c 

J = the instantaneous inertia tensor (deformation dependent) for the body 
The resultant force due to gravity effects is 

(f.) , = 5 * 




( 11 - 118 ) 
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and the forcei acting in the k'* deformation coordinate, is 




adV 


( 11 - 119 ) 


Now, the unit vector, ej^ , has projections onto the body axis system that continually vary 
as the body changes attitude. The unit vector, , is expressed in terms of direction cosines, 
and the three unit vectors associated with the body-reference frame are expressed as 


®R ■ [®bJ I'l'gJ 

Also defined are 

S = [e J {s} 

[sj = SK* 


( 11 - 120 ) 


( 11 - 121 ) 


( 11 - 122 ) 


( 11 - 123 ) 


( 11 - 124 ) 
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With these definitions and the force and torque expressions of equations II-l 17, II-l 18, and 
11-1 19, it follows that the first three elements of the contribution to the right-hand force 
vector due to gravity effects are 

The second three elements are 

{^4 * t H - “’) {*} 

The force due to gravity acting in the deformation mode is 





(11-127) 


where the inertia integrals, N (n = 1, 2, •• 6), and \C^^j (2, m = 1, 2, 3), are con- 

sistent with the development given in paragraph 11. D ana Appendix A. 
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G. Provision for Inclusion of Thermal Environments 

All problems associated with thermally induced deflections have in common the requirement 
that, to determine the effect of solar heating, the spacecraft’s attitude relative to the Sun 
must be known. This required information can be extracted at any point in time from the 
state vector. It is then necessary to have a model of the response of the flexible structure, 
either static or dynamic, to solar heating. 

Considerable work has been done on modeling flexible appendages in thermal environments 
(References 2 through 4), and the results indicate that the response depends on the radiation 
properties of the booms and the attitude relative to the Sun. 

The simulation program accounts for time-dependent thermal deformations in the following 
manner. It is assumed that a model exists whereby the structural deformation of a flexible 
boom (or appendage) resulting from solar heating can be determined from elements of the 
state vector and time. This deformation is subtracted from the actual deformation, and 
the difference is premultiplied by the appendage stiffness matrix. The result is a vector of 
modified, generalized restoring forces for the appendage, which is summed into the 
vector for the appendage body. 

In terms of the development given in paragraphs II.B and II.D where -[k] is seen to be 
the generalized restoring forces (in the deformation coordinates), note that -fk] [ is re- 
placed with-[k] (|^| - thermal deformation state, , is that which must be 

established from a thermal deformation model. 

In this way, a closed-loop response analysis can be achieved, using external subroutines 
to develop the thermal deformations. Some problems may require only open-loop opera- 
tion if the variations of in time is slow with respect to general dynamic response. 

Rather than building a rigid (or irrevocable) model of thermal deformation, the dynamic 
simulation program provides the user with an interface whereby he can formulate and code 
a particular model. Thus, latitude with respect to user requirements is retained. 

III. SYNTHESIS AND ANALYSIS OF THE LINEARIZED SYSTEM 

Developments to this point have described the analytical techniques used to synthesize the 
nonlinear characteristics of a dynamical system consisting of an assembly of interconnected 
flexible (or rigid) bodies. Particular emphasis has been placed on spacecraft systems in 
which individual bodies that comprise the system may be either spinning or nonspinning 
and may have large excursions with respect to each other. 

This section presents a comprehensive summary of the techniques developed for synthesis 
and analysis of the linearized dynamic system with particular emphasis on frequency-domain 
techniques. For the purposes of this discussion, it is convenient to redefine the nature of 
the system under consideration and, in describing the techniques, to consider the total 
dynamic system as a plant subject to a controller rather than as a spacecraft system consisting 
of interconnected bodies that may be subjected to a control system. 
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Linearization of the nonlinear state equations is necessary for applying the powerful analy- 
tical techniques associated with hnear system stability synthesis. When a nonlinear system 
can be reduced to a linear system in the vicinity of a particular state of interest, it is much 
more desirable to work with the linearized state equations. An additional feature related 
to the hnearized system permits the analyst to observe linearized perturbation time-response 
characteristics for the system. The linearized time response can be easily automated by 
recursive formulas, which are generally more efficient than nonlinear numerical integration 
algorithms.* 

A. Introductory Discussion 

The main-line nonlinear time-domain analysis is structured to assemble a collection of inter- 
connected bodies, including a control law. The general form of the governing equations 
may be concisely indicated as 



and the form of the function, F, is the essence of the nonlinear time-domain solution. In 
fact, it can be stated that equation III-l is the fundamental basis for the entire DISCOS 
program. Algorithms for evaluating the nonlinear state-vector time derivatives (and auxiliary 
equations) are centered in a subprogram and its supporting routines. These functional 
algorithms are used for linearizing the governing equations about a specified state. In addi- 
tion, it is desirable to introduce some new variables, including sensor signals, and con- 
trol torques, B. These new variables extend the number of equations, and the additional 
expressions are linearized along with the basic state equations. Additional remarks con- 
cerning the use and manipulation of the additional variables are given in a later section. The 
remainder of this subsection will address specifics relating to the linearization process. 

Attention is first focused on a single variable, y,^, and its dependence on the system state, 

Y‘, through a known (though possibly nonlinear) functional relationship. Arguments begin 
by considering an initial system state, Y’ (o), and a functional algorithm with which to 
evaluate the expression, yj^ = d/dt y,^. The unknown, y^, is first expressed in terms of a 
Taylor’s series expansion about the given state, Y* (o), as 

Yv = Yu (o) + + dYJ dY® + . . . (III-2) 

9Y dY^9Y® 


•Reference Paper I provides a broadbrush narrative description of the theoretical development given in this section. 
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Because the interest here lies in the linear part only, the series is truncated for all partial 
derivatives greater than one and 

dYk 

3YJ 

The task at hand, then, is to establish the partial derivatives indicated as yj^ ^ , thus yielding 
an expression (for all AV = (o), i = 1,2, • • •) of the form 

AY* = H. j AYJ (III-4) 

Because it would be nearly impossible (certainly impractical) to generalize the determination 
of the partial derivatives as explicit analytical expressions involving the independent state 
variables, a numerical approach has been adopted. This task is accomplished by employing 
numerical perturbation techniques in conjunction with quadratic functions to establish the 
desired partial derivatives. Symbolically, determination of elements of H. , is attempted so 
that 


^ = Y‘(o) + HjjAY^ 


(III-5) 


where it is assumed that: 

• The functions, Y’, are indeed linear, sufficiently near the state, Y‘ (o) 

• The functions, Y*, although possibly nonlinear, can be represented as a quadratic 
(or lower order) in the neighborhood of Y* (o) 

The basic approach is concisely summarized in two steps: 

• Establish quadratic coefficients for Y* in the vicinity of the state, Y* (o) 

• Evaluate the partial derivatives, H. ., at the state, Y* (o), using the quadratic co- 
efficients and perturbation values on the independent variables. 

B. Linearization Process 

With reference to figure 8, the quadratic 

f(v) = 


formula can be stated in matrix form as 



(III-6) 
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where 17 is a local spatial coordinate with origin corresponding to and it is desired to 
establish the derivative, 3f/0q, evaluated at q^^^. 

In general, the required partial derivative is 


^ ^ , (in-7) 

dq drj 3q 


Because the three values, ^(i+2)> evaluated by the previously discussed functional 

algorithm, these values satisfy equation III-6. More specifically, consider 



(III-8) 
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and, by matrix manipulation, it follows that 


Kn) 


= Tl 1 j 


— 



-1 



1 

u ) 




1 1 

q i+l 

^i+1 

1 

Vi+ 1 / 

q'i+2 

^i+2 

1 



where the local coordinate, i?, is defined as 


V = 


q - % 
^1+2 ~ 


and it can be noted that 


9t? 1 

= 0: ^i+2 = 1; IT' = 


9q V2 - 


It then foUows that 


f(7?) 


= V ij 


0 0 

q'i+i Vi 
1 1 1 


^-1 
1 

1 


i+l 


i+2 


and, if = 1/2 is specified and it is noted that 


f(7?) - n ij 


— 


— 

y s 

0 

0 

1 

1 ^(0) 1 

1/4: 

1/2 

1 


1 

1 

1 

V(l) / 


(I1I-9) 


( 111 - 10 ) 


(lll-ll) 


( 111 - 12 ) 


(111-13) 
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(IIM4) 



and, in particular, 


E 

dr] 


= f' (t? = 0) = e 

( 0 ) 


and 


3q 


= f'q(4=4|)= 

(°) ^ %^ 2 ) - %) 


(III-15) 


(III-16) 


(III-17) 


1. Comments 

Selection of an initial perturbation value, q(i + 2), from an initial specified state, q(o) = 
Y^(o), is somewhat arbitrary. A value of 1 percent of the initial value has been successfully 
used for all example problems during the course of the study. In the case where the initial 
value is null, an infinitesimal value must be chosen. A value of 1 X 1 0"^ has been accom- 
modated in the digital code. The intermediate choice of T?(i + 1) = 1/2 was selected for 
other reasons. Consider first that a single evaluation of a partial derivative, 8f/3Y‘ is not 
sufficient to qualify its validity. 

An approach has been employed whereby two successive evaluations of 3f/3Y‘ obtained 
by successively cutting the perturbation in half must agree to a predetermined number of 
significant digits (e.g., 5).* The choice of T?(i + 1) = 1/2 r^uires but a single new evaluation 


*The establishment of the ettoi criteria that will be used to compare successive evaluations of 9f/3Y* is strongly dependent 
on computer word size. In special cases, numerical noise can exceed the established error criteria and prevent numerical 
convergence. This problem can usually be circumvented by changing error criteria limits set in subroutine LINEAR. 
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for each element in Sf* at each successive reduction in the perturbation value. In summary, 
the linearization employs an iterative technique to establish the desired partial derivatives. 

2. System Resonance Properties 

The linearization process has provided a system of first-order differential equations that 
describe the dynamical simulation in terms of perturbation variables about an equilibrium 
state. The linearized canonical form appears as 

AYUH,jAYJ (i.j = 1.2,...) (III-18) 


The coefficients, H. . , contain all of the resonance frequency properties of the dynamical 
system. The standard eigensolution form is indicated by taking the transform of this 
expression to obtain 

- Hj j'jAYJ(s) = 0 (III-19) 

Extraction of the roots (eigenvalues) from H. ^ then gives the roots of the dynamical system. 
There will be N of these roots, and any complex roots will appear as conjugate pairs because 
the elements of Hj^. are all real. The imaginary part of the complex pairs represents the 
resonance (or characteristic) frequencies of the system. 

C. Exchange of Variables 

It is often necessary, for the analyst to require additional variables with which to assess the 
stability'characteiistics of the dynamical system. These additional variables ordinarily take 
the form of plant sensor signals and control system output forces and torques. Although 
the desired variables may not be explicitly contained in the system state vector, Y‘, they 
are known in terms of the state variables through an expression of the form 

wJ = g(Y*) (III-20) 

Recall also from previous discussions that it has been established either directly or through 
linearization that 


AY‘=H,jAYJ (III-21) 
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Now, rewriting equation III-20 in matrix form and identifying variables to retain, Yj , and 
variables to eliminate, Y^, gives 



(III-22) 


and it can readily be established that 

hi = w |z| 


where 


[R] = 

and 


Thus, the state equations for the dynamical system can be written (in terms of variables 
that include the desired plant sensor signals and control system forces and torques) as 

|zl = [R]-‘ [HjjJ [r] {z[ (III-24) 

and the transformation A.. = R'^ H. . R, is commonly referred to as a similarity transforma- 

I y 

tion. The matrix. A-::, is said to be the transform of H. . by the matrix R (Reference 5). 

The similarity transformation, A^, possesses a unique property in that the eigenvalues of 
Ay are equal to the eigenvalues of ! A simple proof that establishes this point follows: 

The characteristic matrix of Ay is given by 

(A, - si) = (R-‘ H, , R - si) = R ‘ (H. , - sI) R (III-25) 

tj 1|J 



■C-a 




C. I 




(III-23) 
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It follows that Q(s), the characteristic polynomial of A^, is 

Q(s) = det (A.J - si) = det R'* (det (H. ^ - si)) det R 


and, as (det R'^ ) = 1/det (R), it is apparent that 

Q(s) = det(H. j - si) = P(s) 

where P (s) is the characteristic polynomial of H. .. Thus, it is evident that the matrices, 
Hy and Ay, have the same characteristic equations 

m = P(s) = 0 


and, therefore, the eigenvalues of A., are equal to the eigenvalues of H. .. 

Application of this property now permits isolation of the plant and controller, even for a 
state-space representation of an inherently nonlinear system that can be linearized about 
a specified state. Separation of plant and control-system variables is an important facet 
of linear system stability synthesis. 

1. Evaluation of the Similarity Transformation 

This discussion relates to a procedural approach for determining the similarity transformation 
matrix, [R] , that will relieve the user from the burden of having to select those variables 
to eliminate from the original state vector so that the auxiliary variables, B* and can 
become an independent constituent of the modified state vector for use in the linearized 
studies. With reference to equation III-22, all the Cy coefficients are known because they 
have been obtained through linearization of the auxiliary equations. The Cy coefficients 
simply define the dependence of the auxiliary variables, w>, on the original state variables, 

Y’. In general, it is not possible to directly partition the Cy in the C^ and Cj partitions 
as indicated in equation III-22, for the decision has not been made yet as to which state 
variables to retain and which to discard in preference to introducing the auxiliary variables, 
w*. In this light, a best possible choice is made with regard to which of the variables to 
eliminate from the state vector, Y*, so that the auxiliary variables, w^, may be included. 

A one to one variable exchange will often occur between an element of w^ and an element 
of Y*. In any case, a variable exchange is necessary for structuring the total system into 
the desired plant/controller framework whereby the plant and controller can be isolated 
along with the plant sensor signals and the control-system inputs. 
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The following approach is used in this simulation to accomplish the desired result (namely, 
an optimum selection from Y* as to which variables to eliminate so that w^ can be introduced 
as a part of the state vector). With reference to equation III-22, 



aiI-26) 


The primary focus of attention is now directed to a systematic examination of the co- 
efficients so that the variable exchange is accomplished in an optimum manner. To help 
clarify the discussion, some size identifications are noted: 

Cjj has size NR by NS 

Y* has size NS by 1 

w> has size NR by 1 

NJQ = NS + NR 

Clearly, at least one nonzero element exists in each row of the C.j array. Otherwise, Y* does 
not represent an independent set. 

Now, a search through the first NS elements of row 1 in the matrix array 



(III-27) 


will identify the largest element (absolute value) in row 1 . Assuming that this element occurs 
in column JBIG (1 < JBIG < NS) permits the division of each element of row 1 by this 
largest element, and subsequent elementary row operations on rows 1 through NR will 
eliminate those elements below the pivotal element in column JBIG. This procedure is 
repeated for each of the NR rows contained in the matrix, and the following observations 
are noted: 

• The appearance of a one (1.0) in a row identifies a variable that will be eliminated 
in preference to including an element of w^. 

• The absence of a zero or one in columns of a given row indicates which variables 
will survive the exchange process. 

• All variables in w^ (NR of them) will become part of a new and independent state 
vector (the modified state vector). 

• The transformation, Rjj (i, j = 1 . . . NS) can be constructed from the matrix that 
remains after the procedural approach has exhausted all of the NR rows of expres- 
sion III'27. 
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2. Illustrative Example 


A simple example is presented to further describe the actual mechanics used for evaluating 
the similarity transformation. Linearization of the auxiliary equations has established that 



and, to implement the proposed algorithm, the C matrix augmented with -I is first written 
and the elements of C are further identified as 









1 ^ = 

^11 ^2 ^13 *^14 

O 

1 

1 — 1 
1 

*^21 ^22 ^^23 *^24 

\ 

1 

o 


For illustrative purposes, assume that I bjj l>l (bjj, b^j, bj^) I, and, hence, b^j is the 
pivotal element for the first row. It follows that row 2 is modified by the relation 


= - ^23 


+ b. 


13 


giving the matrix 


\ 




1 

bl4 

1 

n 

b,3 

bl3 

^3 

"^3 

U 

b.i 

K — — 4. K 

bl2 

-h - j. K 

0 

bl4 

h h h 

^23 

-1 


°23 u ^ *^22 

'’l3 

” b,3 

bi3 


\i 

S2 

1 

S4 

Ss 

o' 

Sr 

S 2 

0 

S4 

Ss 

-1 
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The process is continued by first dividing row 2 by the largest element and again using row 
operations to eliminate the other element in that column. Again, for illustrative purposes, 
it is assumed that the pivotal element of row 2 (say C 22 ) is known, and row 1 will be modi- 
fied by row operations as 


ele,. 


12 


^2j 


22 


+ C, 


giving 


Si 

0 

1 

*^24 

Ss 

S 2 

■S 2 ~ *^11 

S 2 

1 

^22 

“S2~ Ss 

^22 

S 2 

Si 

1 

0 

S 4 

Ss 

1 

S 2 

S 2 

*^22 

1 

[ 

r ^ 

1 


dn 

0 

1 

dl4 

dis 

dl6 

^21 

1 

0 

^24 

^25 

^26 


from which the desired similarity transformation is established as 


1 

0 

0 

0 

-^i21 

-^24 

-‘ 12 S 

-^26 

-dn 


-dis 

-di6 

0 

1 

0 

0 
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It, now follows that the original state variables are written in terms of the modified state 
variables as 



D. System Transfer Functions 

The entire system transfer function synthesis can be concisely summarized in a chronological 
sequence of steps that began with linearization of the coupled mechanical/control law 
equations that govern the dynamical motion. This process included linearization of addi- 
tional equations that contained specific variables required for further consideration in the 
stability analysis (namely, plant sensor signals and control system outputs). A similarity 
transformation has been introduced that exchanges original state variables for these desired 
sensor signals and controller outputs so that the resulting modified state vector is still repre- 
sentative of an independent set of state variables. The resulting system of state-space 
equations is later identified as equation III-28. 

The system characteristic matrix, Ajj, provides the basis for evaluating the coupled mechanical/ 
control system resonant characteristics (natural frequencies), as well as the fundamental basis 
for specification and determination of the various types of transfer functions. The next sub- 
section addresses some of the more specific details regarding specific transfer function relation- 
ships. A particular transfer function is identified by a type along with the desired output/ 
input variable designation. An eigenvalue problem is then stated, which leads to determi- 
nation of the numerator roots (zeros) and denominator roots (poles) for the particular trans- 
fer function. When the poles and zeros are known for a transfer function, this information 
can be further processed and displayed by any of the conventional display modes: Bode, 
Nichols, Nyquist, and/or root locus. 

The conventional block diagram representation for the coupled plant/cbntroller system 
(figure 9) provides additional insight for determination of system transfer functions. 

The first-order differential equations for the system are written as 


V = A., ZJ + B_ R_J + B R J 

ij T,, T Sj, s 


(III-28) 
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Figure 9. Plant/controller block diagram. 


and it is helpful at this point to express the equation in matrix form and indicate the sepa- 
rate partitioned subsets of Z‘, Ay, Z^, and as 



The following observations are noted 


Si 

= 0 

b.j.1 

= -S4 


= 0 

Si 

= 0 

b.j.2 

^24 


= 0 

^13 

= 0 

b.j.3 

= 0 

b 3 
s 

^32 

Sa 

= 0 

bx4 

= 0 

b 4 
s 

cd 

It 
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and equation III-29 can be restated as 




(III-30) 


Equation III-30 is the operating basis for stating particular transfer function relationships 
for the plant/controller system. 

The general procedure is to establish a system transfer function between inputs Rj and 
and outputs and B. Loops may be opened to provide open-loop information by mani- 
pulating the Ay coefficients to prohibit certain feedbacks. 

To symbolically describe specification of a transfer function, begin by consolidating the b 
coefficients and taking the Laplace transform of equation III-30 to give 


or 



(III-31) 




(III-32) 


and then employ Cramer’s rule to evaluate a given element, Z(sF , due to a particular input, 
U(s)’, where 


O'MS) 

and where aug I Is - A I is accomplished by placing column q of b into column p of I Is - A I 

The Q-R algorithm (References 6 and 7) is a useful tool with which to extract the indicated 
determinants in equation III-33. 
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1. Root Extraction Process 

With reference to equation III-33, evaluation of both the numerator and denominator roots 
is desired. The denominator root extraction is straightforward in that , P 2 > P 3 > * ** Pn 
is sought from an expression of the form 

D(s) = det(p]s - [A]) 

so that 


D(s) = (s - Pj)(s - Pj) • • • (s - P„) = TT - Pj) (III-34) 

i=l 


This evaluation is completed by extracting the characteristic roots of the matrix, Ajj. In 
general, these roots will be complex because Ay is not symmetric. 

The process used for evaluating the numerator is best illustrated with an example. Consider 
that we have the (4 by 4) characteristic system matrix, 


^11 

S 2 

S3 

S4 

Si 

^22 

^23 

^24 

Si 

^32 

^33 

^34 

Si 

S 2 

S 3 

^44 


and the column of coefficients bj that premultiply the desired input variable, . Further, 
let it be desired to obtain the transfer function relating output of the third variable in the 
state equations, yj , to the input, U^. 

The state equations for this system would appear as 



(III-35) 


68 



and, with reference to equation III-33, the numerator is 


N(s) = aug I Is - A I 


or 


N(s) = det 


s-a,, 

-S2 


-“l4 




-S4 

-Si 

-®32 


-S4 

-Si 

~^42 


S-344 


(III-36) 


After performing elementary row operations, equation III-36 can be restated in the form 



* -Si 

Si 

-Sa + Sa^’i/^ 

-S4 + S 4 ’^l/'’3 

N(s) = bjdet 

-Si + 

Si ‘’2/‘^3 

S-S 2 + S2^/’>3 

-S 4 + ^34 


-Si + 

Si ‘’ 4/^3 

-S 2 S 2 *’ 4/^3 

S-344 + S 4 ’’4/*’3 


(III-37) 


or, in symbolic terms, as 


N(s) = bj det 1 [Is] - [a]i| 


(III-38) 


where the matrix a is given as 


Si -Si 

^12 "Sa ^1^^3 

S4 -S4 ’’ 1/^3 

Si -Si 

Sa-^saV^ 

^24 -^34 '^2^^3 

Si ~^31 S^^3 

S 2 -S 2 N/^3 

“44 -“34 S/*’3 


Note that the previous 'expression for N (s) is finite only if b 3 # 0, and the question is : Can 
bj realistically be null? The answer is yes, as the following example indicates. 
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Consider the simple mechanical system consisting of two masses connected by a single spring/ 
dashpot combination as shown in figure 10. 



Figure 10. Simple mechanical system. 


The state-space representation is 



and the frequency domain (or transformed) equations in s are 
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where 


[A] = 

-c/nij 

c/nij 

-k/mj 

k/nij 


c/m^ 

-c/nij 

k/nij 

-k/m. 


1 

0 

0 

0 


0 

1 

0 

0 


Now, consider the transfer function, Xj(s)/Fj, where the augmented numerator is 


N(s) = det 


l/nij 

-c/mj 

k/m^ 

1 

3 

0 

s + c/nij 

-k/m^ 

k/mj 

0 

0 

s 

0 

0 

-1 

0 

s 


and the pivot element is the (1, 1) element or 1/m^ ^ 0. On the other hand, the transfer 
function, Xj (s)/F^ , has the augmented numerator 


N(s) = det 


s + c/mj -c/mj 1/nij -k/nij 
-c/m^ s + c/mj 0 k/nij 
-10 0 0 

0 -10 s 


and the pivot element is the (3, 3) element, which is null. 

The problem to be addressed nowiinvolves evaluation of the numerator determinant, N(s), 
when the pivotal element is null. The particular mathematical problem may be restated as 


N(s) = det 


[T]s - [A], 


(III-39) 
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where [I ] is the identity matrix, [I] , of size N with a null diagonal element. Addition and 
subtraction of the quantity, [I ] x > where x is an arbitrary constant not equal to one of the 

fsf 

roots of [A] , yields 


N(s) = det 


[I] (s - x) - [A] - [I] X 


[[A] - [T] x] 


(III- 40) 


and, if (s - x) = 1/p is defined, there results 


(-1^ 

N(s) = det 


/>W 


l( 


A - Ix det Ip-(A-Ix) ‘l 


I) 


(III-41) 


The roots, (Pj, i = 1, N), are found as the eigenvalues of the expression 


[[A] - [T] x] [T] 


(III-42) 


and the eigensolution permits N(s) to be written as 


(-l)N 

N(s) = det 


A - Ix 


(p - Pi)(p - Pj) • * * (p - Pn)> (in-43) 


The following observation can now be made: a Pj equal to zero implies a root at infinity 
(or a characteristic polynomial having an order less than N). Thus, the null Pj’s are elimi- 
Inated from the expression, giving the characteristic polynomial an order, n, which is less than 
N. It is a rather common occurrence for the number of zeros (an order of N(s) to be signi- 
ficantly less than the number of poles (order of D(s)). With reference to equation III-43, 
the numerator expression, N(s), can be written as 


N(s) = (-1)^ det 


A - Ix 




(III-44) 
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and, recalling that p = 1/s-x, yields 


i-lf 

N(s) = — det 

fr - V 

i=l 


A - 7x |(s - Sj)(s - Sj) • • •■(s - 


(III-45) 


Next, note that 


1=1 


TT (x-s,)=TT 

i=i \ Pi 


(III-46) 


and it follows that 


N(s) = (-1)^'" Tf pj det ; A - Ix ; TT (s - s.) 


i=l 


i=l 


(III-47) 


The numerator root gain, kj^ , can now be identified as 


kR =(-!)''■" TT Pi det 


1=1 


A - Ix 


(iII-48) 


and the Bode gain, , for the numerator is 


kg = kg (-O'" TT 


i=l 


where m < n 


2. Transfer-Function Classification 

With reference to figure 9, it is possible to directly identify six transfer function types, each 
characterizied by the specific variables involved and by the presence of feedback. In addition, 
a seventh type will also be described whereby certain control variables feed back and others 
do not. This type is similar to an open-loop transfer function but treats selected channels 
of the controller as part of the mechanical system (plant). During the course of this discussion, 
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it will become apparent that additional transfer function types are easily accommodated by 
rather simple manipulations with the system characteristic matrix, A^.* 

In general, note that the process of obtaining the desired transfer function involves but a 
few basic steps. The transfer function characteristic matrix, IRjj, and the desired force 
coefficient vector, bj, are obtained directly from the system characteristic matrix. Ay, These 
two matrices are then put in a form so that the Q-R algorithm can be used to extract system 
roots. 

a. Type / (Plant On/W — Type I is the forward path transfer function for the plant with no 
feedback and is of the form 


= G(s) (III-49) 

The control variables, 5‘, and control outputs, B‘, do not feed back into the plant. The matrix 
expression depicting the system of interest is 


£ 

dt 



(III-50) 


The matrix. Ay, to use in the general expression given as equation III'33 is referred to as 
fRy or the reduced Ay matrix, 


fR. 




Si S2 


(III-51) 


The augmented fRy matrix is obtained by removing the column corresponding to the input 
variable, R^, from the expression, b.j., and inserting this column into the column in fRy, 
which corresponds to the desired output, The resulting transfer function is then given as 


X„'/R’ = 


aug I Is - fR| 
I Is - IRl 


(III-52) 


‘Additional comments pertaining to the exact definitions of various types of transfer functions implemented can be found 
in subroutine DEF5 Debug- 1 1 6 in the discussion concerning “Frequency Domain Analysis ” (Volume U, Appendix B, 
page B-332). 
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b. Type // (Controller Qn/W — Type II represents the feedback path, H(s), for the controller 
only. The desired transfer function relates control-system outputs, B', to sensor signal in- 


BP/X^ = H(s) 


(III-53) 


The reduced characteristic matrix, fRy, and the corresponding input coefficients, b.^, are 
given as 


II 

^33 

1 

> i>ik = 

^32 


^43 

^AA 

44 


h2 



• 


_ 


(111-54) 


c. Type III (Open Loop, /yp — Type III falls within the framework of the classical open-loop 
transfer function designation and relates control system outputs, B*, to external plant inputs. 
Rip. The algebraic expression for a given output variable, BP, due to an external input, R^, 
is indicated as 


BP/R^ = (HG)(s) (UI-55) 

The open-loop system characteristic matrix, fRy, and corresponding input coefficients, bj^, 
are 



^12 



It 

1 

^21 

^22 




-“24, 


^32 


“43 


0 


^42 

^43 

“44 


0 


(III-56) 


It was previously noted that ajp = a^p ^ 2 ^ ~ ^ addition, the partitions, aj^ 

and aj 4 , are set to zero to prohibit the B* feedback. Thus, the loop is opened to establish 
HG, and the open-loop transfer function in s. Note that the negative sign in the bjj^ coeffi- 
cients simply indicates that the B' feedback is negative with respect to the external plant 
inputs. Rip. 
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d. Type IV (Open Loop, GH j—An additional open-loop transfer function is often desired 
for assessing the plant sensor signal outputs caused by controller noise inputs. The transfer 
function then relates sensor signal outputs, X’ to control system noise inputs, . The 
plant sensor signal vector does not feed back into the system so that 

X = (GHXs) (III-57) 

and the system characteristic matrix, rRy, and the external input coefficients, b.^, are identi- 
fied as 


Si 

^12 


S4 

II 

1 " 
o 

\ 

Si 

^22 


®24 


0 




^a4 


S2 




^44 


S2 


(111-58) 


Note that the a ^2 a ^2 partitions have been nulled to eliminate sensor signal feedback. 

e. Type V (Closed Loop — Control Ratio) — The system control ratio is given as the transfer 
function that relates plant variable outputs to externally applied plant inputs with the control 
system entirely active. This transfer function is expressed as 


Xss‘’/n = 



(III-59) 


and the system characteristic matrix, fRy, and the external input coefficients, bjj^, are identi- 
fied as 


Si 

S2 


S 4 

’ '’ik = 

- — J 

Si 

S2 


S4 


-S4 


S2 

Sa 

*a 4 


0 


S2 

Sa 

®44 


0 

_ 



__ 




(111-60) 


The negative sign in the matrix, bjj^ , indicates that the feedback is negative. 
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f. Type VI (Closed Loop) —kn additional closed-loop transfer function has been accommo- 
dated within the digital simulation. Specifically, Type VI relates plant sensor signal outputs 
to sensor signal noise inputs with all control-system loops active. The transfer function is 
symbolically indicated as 


= (transfer function) (III-6 1 ) 

where the system characteristic matrix, IR^, and corresponding input coefficients are identi- 
fied as 


^11 

S2 


1 


( 

0 

i 

Si 

S2 


S4 


0 


S2 

^33 

S4 


®32 


S2 

S3 

^44 


S2 






• • 


(III-62) 


g. Type VII ( Quasi -Open Loop) -An additional transfer function type is identified here and 
is referred to as a quasi-open loop (figure 1 1). It is of the open-loop type in that we are 
interested in control-system outputs, B‘, attributable to plant variable inputs. Rip. For 
example, suppose that, for a multichannel control system (such as azimuth and elevation), 
outputs B* are desired on the controller channel that do not feed back and suppose that the 
other channel is active in that it feeds back into the plant. 



Figure 11. Type VII quasi-open loop block diagram. 
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For the configuration indicated, a typical Type VII transfer function (TF) would be given by 


B| = (transfer function) 


and the form of the system characteristic matrix, IRy, and plant input coefficient matrix, 
bj^, would be 


fR ij 


'‘11 

S2 


S4 


-*14 

Si 

S2 


^24 


-*24 


S2 

^33 

*34 


0 


S2 

S3 

*44 


0 


(III-63) 


The subpartitions, a^^ and 3 .^ 4 ^ indicate modification of the original partitions, aj^ and 324 - 
Specifically, a^^^ is a subset of a., that is obtained by keeping only those n columms of aj^^^ 
that correspond to the B' variables that feed back to the plant. 


3. Transfer Functions— Polynomial Description 

This subsection is addressed to implementation of control-system transfer functions described 
as the ratio of two polynomials in the frequency domain, s. Specifically, consider 


TF = P(s)/Q(s) 


(III-64) 


where 


Q(s) = + 3jS + + • • • + 


and 


P(s) = bjj + bjS + + 


m 
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Because the previously described governing equations have been stated in canonical first- 
order form, the polynomial description for the transfer function is restated in the form 


S‘ = + B,U 


(III-65) 


The block diagram for the system is 


u 

P(s) 

s 


Q(S) 



from which we write 


P(s) 

6 = — ^ U 

Q(s) 


( 111 - 66 ) 


and expansion of the implied operator in s results in a differential equation of the form 

6 ‘ + . . . + I + U + b^.j u'* + • • • + b, U + b^U (HI-67) 

where 


5 - — 

dt" 


In general, the order of P(s) will be no greater than the order of Q(s) or m < n. 
a. m =77— Equation III-67 is divided by a^^ to obtain 


n n-1 • m m-1 

6 + C , 5 + + C, 6 + C„6 = d„ U + , U + 

n-1 10 m m-1 


+ djU + d„U (III-68) 


where 


C| = 
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and 


b. 

1 



The following example is used for illustration: 
Consider the equation with m = n = 4, 


5 + C3 6 + CjS + CjS + Cq 5 = U + dj U + d^U + djU + d^U 


or, in operator form, 


+ s^CjS + s^CjS + sCj6 + CqS = s'^d^U + s^dgU + s^d^U + sdjU + d^U 


This can be rewritten as 


(5 - d,u) + s3(c 38 - d3u) + s2(c35 -djU^) + s(c,6 - d,u) + (c^S - d„U ) 


and the substitution 

5j = 5 - d^U 

permits a reduction in order to 

8^(5^ + C35 - d3u) + s^{c ^8 - d^u) + s(c ,6 - d,u) + (C58 - d„u) = 0 

A new variable can again be introduced: 

5, . (5, +C,6 -d,u) 
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and the foregoing can be rewritten as 


s" (6^ + C^a - d^U) + s(Cj5 - djU) + (C„6 - d„U) = 0 


It follows that, if 

«3 = S*2 - d^U 

is defined, 


s( 53 + Cj6 - d,U) + C„6,- d„U = 0, 


results, and the substitution 


a 


4 


= h ^ - djU 


gives 


«4 = -Co« 


+ d„U 


The variable, 5, can now be eliminated from each of the foregoing expressions, and the re- 
sults can be generalized to n'* -order systems. 

The result is concisely stated as a matrix equation that is recognized to be of the desired 
form initially given as equation III-65, 
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where 6^ and S,!the original variables of the equation, are related as shown previously, and 
U is the input variable to the transfer function expression as indicated in equation III-65. 

b. m < n— The general expression for the case where m < n is easily accommodated by 
restricting the dj coefficients to reflect the limit m. Commonly, only the d^ coefficient 
will be finite. 

4. Frequency Response 

Transfer-function poles, zeros, and root gain can be converted to the standard Bode form 
for frequency response by combining time constants, damping, and resonant frequencies as 


TF = 


N1 


N2 


2f,S 


TT (1 + Tjs) TT (i + — + — 


i=l 


Ml 


i=l 


J=1 


M2 


TT (1 + TT 


CO. 

1 


j=l 1 + 


CO. 


CO. 


2fjS 
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where the Bode gain is 


where 


kg = k 


TT 

i=l 


Z. 


m 

TT P) 

j=l 


k = root gain 
T = system constants 
f = system damping at frequency co 
CO = system resonant frequency 

The frequency response is then calculated by substituting jco for s and evaluating the transfer 
function expression at various co’s. The digital simulation uses a vernier frequency incre- 
menting approach that automatically introduces smaller frequency increments near the poles 
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and zeros. This variabie frequency incrementing technique permits better transfer-function 
resolution near the resonances at which amplitude and phase can vary rapidly.* 

5. Root Locus 

The root-locus method of analysis and design is based on the relationship between the poles 
and zeros of the closed-loop transfer function and those of the open-loop transfer function. 
The method is used to determine the location of the roots of the characteristic equation 
as a function of a single open-loop gain parameter. The locations of these roots are indicative 
of the relative system stability. The analyst may use the method as a design tool by adjust- 
ing the poles, zeros, and the open-loop gain parameters so as to yield a closed-loop system 
with satisfactory critical frequencies (poles and zeros). 

To further describe the theoretical basis for the method, refer to the conventional control 
ratio for a feedback system as shown in figure 12. 



Figure 12. Conventional feedback control system. 
The control ratio, C(s)/R(s), is 


C(s) ^ G(s) 

R(s) l + G(s)H(s) 


(III-71) 


♦Eigenvector analysis in coi\junction with frequency domain analysis is not a common practice; however, for complex, 
highly coupled multirigid and flexible-body systems, it becomes a necessity. System roots often shift significantly from 
body and control &equencies, and are therefore impossible to physically interpret. Eigenvector analysis of the character- 
istic matrix associated with the transfer function under study can be used to determine a measure of the degree to which 
each of the body and control-system degrees of freedom couple to form the total system roots and associated modes 
(eigenvectors). For further information, see subroutine DEF5. 
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and the open-loop transfer function, G(s) H(s), is identified as a ratio of two functions in s, 


G (s) H (s) = k 


P(s) 

Q(s) 


(III-72) 


The characteristic system equation is 

1 + G (s) H (s) = 0 


(111-73) 


1 + k 


P(s) 

Q(s) 


= 0 


The conventional root-locus plot portrays the loci of the values of s that satisfy the character- 
istic equation as k varies from zero to infinity. Note that 

• At k = 0, the roots of the characteristic equation are equal to the roots of Q(s), 
which are the same as the poles of the open-loop transfer function, k P(s)/Q(s). 

• As k approaches infinity, the roots approach the roots of P(s), the open-loop zeros. 

Thus, as k varies from 0 to infinity, the loci of the closed-loop poles migrate from the open- 
loop poles to the open-loop zeros, and the direction of migration depends on the sign of 
the open-loop gain parameter, k. 

Rewriting equation 111-73 yields a more conventional expression for the characteristic 
equation as 


k 

Q(s) 


= -1 


and two conditions are required as follows: 


(111-74) 


Q(s) 


/P(s)/Q(s) = 180°, k> 0 
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The first of these conditions can be expressed as 


Q(s) 

P(s) 


for those values of s that satisfy the angle criterion. The conditions that govern the migration 
of the roots in the complex plane can be solved by an iterative procedure. The iterative pro- 
cedure for evaluating* a single root locus* is described in Appendix E. 

E. Linear Time-Domain Response 

The linearized canonical first-order system of equations can also provide a basis for studying 
system time history in terms of perturbations about a specified state when the system indeed 
behaves in a linear manner in the vicinity of the state. The nonhomogeneous form of the 
equations, the basis for determining system transfer functions, appeared previously as 

Z‘ = AyZl + bj^U’^O) (III-75) 


The external system inputs are the elements of It is convenient to establish the solution 
for the foregoing system through the use of a recursive-formula numerical-integration proce- 
dure rather than throu^ the Runge-Kutta approach. 

Consider the Adams’ corrector formula (Reference 8) at time, t + 1 , 


't+i 


TJ + 

* 24 


9Vi + 19 




(III-76) 


where h is the incremented time step. 

Application of this formula to our system of equations gives 


z* — 7 } 
^t+1 ^ 


24 


j^9Ayz{^, + \9’z\ - 5zj_j + z-.^H 


‘Welch, Raymond V. NASA/Goddaid Space Flight Center, Branch Report No. 254, October 2, 1973. 
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and manipulation yields the solution for all the z* at time step, t + 1 , 




Note the requirement for z* at timeistep, t - 2; hence, the requirement for a starter (e.g., 
Runge-Kutta) for initiating the solution process. 


Goddard Space Flight Center 

National Aeronautics and Space Administration 
Greenbelt, Maryland October 1977 


86 



REFERENCES 


1. Likens, P. W., “Analytical Dynamics and Nonrigid Spacecraft Simulation,” Technical 
Report 32-1593, Jet Propulsion Laboratory, Pasadena, California, July 15, 1974, 

2. Fixler, S. Z., “Effects of Solar Environment and Aerodynamic Drag on Structural 
Booms in Space,” J. Spacecraft, 4|(3), March 1967. 

3. Fiisdi, H. P., Coupled Thermally Induced Transverse Plus Torsional Vibrations of a 
Thin-Walled Cylinder of Open Sections, NASA TR R-333, March 1970. 

4. Goldman, R. L., “Influence of Thermal Distortion on the Anomalous Behavior of a 
Gravity Gradient Satellite,” AIAA Mechanics and Control of Flight Conference, 
Paper No. 74-922, Anaheim, California, August 1974. 

5. Hovanessian, S., and L. A. Pipes, Digital Computer Methods in Engineering, McGraw- 
HUl Book Company, New York, 1969. 

6. Francis, J. G. F,, “The QR-Transformation— A Unitary Analogue to the LR-Transfor- 
mation,” The Computer Journal, 4, Part 1, October 1961. 

7. Francis, J. G. F,, “The QR-Transformation-A Unitary Analogue to the LR-Transfor- 
mation,” The Computer Journal, 5, Part 2, January 1962. 

8. Schied, F., “Theory and Problems of Numerical Analysis.” Schawn’s Outline Series, 
McGraw-Hill Book Company, New York, 1968. 


87 



REFERENCE PAPER I 


DYNAMIC RESPONSE AND STABILITY ANALYSIS 
OF FLEXIBLE MULTIBODY SYSTEMS 
(from DYNAMICS AND CONTROL OF 
LARGE FLEXIBLE SPACECRAFT, 
Proceedings of a Symposium held 
in Blacksburg, Virginia, 

June 13-15, 1977) 


Preceding Page Blank 



DYNAMIC RESPONSE AND STABILITY ANALYSIS OF 
FLEXIBLE, MULTIBODY SYSTEMS 

by 

Carl S. Bodley and A. Colton Park 
Martin Marietta Corporation, Denver, Colorado 
A. Darrell Devers 

Texas A and M University, College Station, Texas 
Harold P. Frisch 

Goddard Space Flight Center, Greenbelt, Maryland 


Abstract 


An approach for d 3 mamlc sinoilatlon and stability analysis of systems 
of interconnected flexible bodies is discussed. The overall approach is 
unique in that any member body of the system can be flexible and the 
total system is not restricted to a topological tree configuration. The 
equations of motion are developed using the most general form of La- 
grange's equations including auxiliary nonholonomic, rheonomic conditions 
of constraint, Lagrange multipliers are used as interaction forces/ 
torques to maintain prescribed constraints. Nonlinear flexible/rigid 
dynamic coupling effects are accounted for in unabridged fashion for in- 
dividual bodies and for the total system. Elastic deformation can be 
represented by normal vibration modes or by any adequate series of 
Rayleigh functions, including so-called quasi- static displacement func- 
tions. 


A digital computer program system has been developed to numerically 
implement the modeling and analysis capability for a wide range of dyna- 
mic simulations including the nonlinear time domain and/or linear fre- 
quency domain. In particular, application has been made to: (1) time 

domain solution for nonlinear response of systems idealized as a collec- 
tion of individual bodies; (2) numerical linearization of system govern- 
ing equations; (3) time domain solution for perturbation response about 
a nominal state; and (4) frequency domain stability analysis correspond- 
ing to the linearized form. 


Introduction 

State of the art dynamic response analysis of a system of inter- 
connected bodies, typical of current and projected large spacecraft, has 
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been restricted until recent years, to rigid bodies interconnected in a 
topological tree configuration. The formalisms of Hooker and Margulies 
(Reference 1) and Roberson and Wittenburg .(Reference 2) set the prece- 
dent with respect to a procedural synthesis and solution of the equations 
of motion governing such a system of bodies. Conq>uter programs, with 
associated formalisms, have also been developed by Velman (Reference 3), 
Frisch (Reference 4) and Farrell, Newton and Connelly (Reference 5). 

These have proven to be highly applicable and useful tools for studying 
spacecraft dynamic response in the time domain. Following these earlier 
formalisms, a host of investigators (for example. References 6, 7, 8 and 
9) extended the earlier concepts to include flexible terminal bodies and 
the capability to deal with environmental effects such as gravity gradi- 
ent torques. 

Almost all current spacecrafts are configured such that existing 
analytical techniques are conpletely adequate. That is, a t3q>ical space- 
craft can be modeled as a rigid central body with flexible texminal bodies 
interconnected in a topological tree configuration. However, it is anti- 
cipated that future spacecraft designs will include configurations that 
cannot be adequately represented by a topological tree. To be more 
definitive, a non- topological tree configuration is one wherein there 
exists re-entrant branches or closed loops of connectivity. Such a sys- 
tem has, in effect, a statically indeterminant load path; this creates 
special problems for a dynamic analysis formalism. 

This paper describes, in summary form, an analytical technique and 
associated computer program system which deals effectively with the 
problems of re-entrant connections. The analytical considerations and 
DISCOS* computer program system are documented in detail in Reference 10, 


Definition of the Physical System 

The spacecraft, which is subjected to analysis, can be described as 
a cluster of contiguous, flexible structures, or bodies. In the analy- 
tical treatment, the entire spacecraft system, or portions thereof, may 
be spinning or nonspinning. Member bodies of the system are capable of 
undergoing large relative excursions such as those of appendage deploy- 
ment, or rotor/stator motions. Bodies of the system may be interconnected 
by linear or nonlinear springs and dashpots; they may be interconnected 
via a mechanism that consists of gimbal and slider block or a servo- 
actuator, or any combination of the above. 


Analytical Considerations - Nonlinear State Equations 

In conjunction with this discussion that follows, the Interested 
reader is urged to see Bodley and Park (Reference 11) where attention 


* D3mamic interaction and Simulation of Controls and Structure 



is concentrated on development of governing equations of motion for a 
single flexible body. In the referenced development, ordinary momenta 
components (as opposed to the generalized momenta of Hamilton's canonical 
equations) are used as state variables. The ordinary momenta correspond 
to nonholonomic velocities (or, time derivatives of quasi-coordinates) 
in the same way that generalized momenta correspond to generalized velo- 
cities. The development has been used by the authors to study dynamic 
behavior of a single flexible body and also it appears to have withstood 
the test of time with respect to general validity and accuracy. 

Presuming the validity of that development, the equations of dynamic 
equilibrium are 


{p} = {g| (1) 

where |p| represents a column vector of ordinary momenta and where |g[ 
represents a collection of all external forces, internal restoring and 
dissipative forces, state dependent inertial forces (i.e., gyroscopic, 
Coriolis, Euler forces, etc.), and servo- actuator forces/torques. Devel- 
opment of the explicit form of the |g} vector is very lengthy and is, 
therefore, not given here. The interested reader will see Reference 10 
where the various contributions to the |g| vector are treated with indi- 
vidual detail. Let it suffice to comment that it is possible to account 
for all state dependent and explicitly time dependent effects in an un- 
abridged fashion. 

Now, the vector of ordinary momenta, |pj is related to a vector of 
nonholonomic velocities |u| through a deformation dependent mass matrix 
[m] ; that is. 


{p} = [™] {l>} , (2 

Ti^iere the elements of [m^ generally depend on the displacement coordin- 
ates |||. Thus, the mass matrix has a time derivative. Implementation 
of the equations of motion of Reference 11 for solution via digital com- 
puter did not require evaluation of [m] and this was one attractive fea- 
ture of the formalism. Such is not the case when considering a cluster 
of interconnected bodies, when the formalism uses Lagrange multipliers 
(interconnecting constraint forces/torques) and kinematic equations of 
constraint, as will be seen. 

The equations of dynamic equilibrium for the j th body of the system 
may be extended to include, on the right hand side, the effects of the 
interconnection forces and torques and thus, 

■ Hi * [■’JII'I- 
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Also, the equations of kinematic constraint are easily expressed as 


E[‘l, W.-I'l 


where [^] j represents a matrix of kinematic coefficients for the j th 
body and \mere |a(t)} is a column of prescribed time dependent velocities 
across the interconnection junctions. The subscript j ranges from 1 
through the number of bodies of the cluster. 

Now, to be able to implement Equations (3) and (4) in a numerical 
integration solution process, it is necessary to numerically evaluate the 
highest derivative elements; in this case, evaluate |p[ , given initial 
values for |p| . However, before this is possible, it is necessary to 
view Equations (3) and (4) as simultaneous equations involving |p[ and 
j X I as sets of unknowns. To simultaneously satisfy Equations (3) and 
(4), let us differentiate Equation (4) with respect to time and change 
the variables of Equation (3) as follows: 

{p} - -s (["]{"})• 

thus, 

l“}j ■ ["]j 

H- EH. H.- 

j j 


(5) 


(6) 


Now, substituting Equation (6) into Equation (7) yields 

H ■(£[“]. [-L'M'nr 

•[•].[■];■ MJ 

L j 


( 8 ) 
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The numerically evaluated | X | given by Equation (8), when substitu- 
ted into Equation (6), provides complete numerical definition of |u| , 
v^ich represents a vector of time derivatives of nonholonomic velocities 
(accelerations as seen by the observer moving and rotating with the ref- 
erence frame of the j th body). The collection of |u|j vectors for the n 
bodies of the cluster serves as a portion of the complete state vector. 
Now, as the, expression to evaluate the | X } and | G* | vectors and the [m] j 
matrix must involve dependence on the elements of {u [ j and on position 
and displacement coordinates | /3 | and | ^ [j » there must be additional 
differential equations to provide | )3 J and | feedback. These are of 
the form; 


j 

{«} J ' [^]j {“}j 


where | 0 | represents a vector of position coordinates (gimbal angles and 
cartesian position vector components) and where | f}j represents dis- 
placement coordinates corresponding to assumed displacement functions 
(eigenfunctions, or any other kinematically admissible functions that 
are consistent with the boundary conditions). In Equation (9) the [b] j 
represents a matrix of kinematic coefficients of very similar form to 
the [b]j of Equation (6). That is, the elements of [b ] j and [ B ] j are 
each dependent on the elements of | 0 } and | f | . Finally, in Equation 
(10), the [s]j matrix simply provides a selection operation as | | is 
a subvector of ]u [j. 

The collection of | vectors, for the n bodies, and the | /3 | vec- 

tor are truly the generalized coordinates of the structural system (they 
provide a complete configuration description of the system of bodies) . 

We note that the equations of dynamic equilibrium (6) involve generalized 
coordinates and velocities in an implicit fashion, while generalized ac- 
celerations are not evident. 

To complete the definition of the state equations, one must account 
for additional differential equations such as those linking control 
actuator torques to combinations of state elements (attitudes, their 
rates and time integrals) through an appropriate control law. These are 
of the form 

{ 5 }= f { U|, t). (11) 
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The explicit form of the "controller". Equation (11) , must be pro- 
vided by a prospective DISCOS program user in that there is an unlimited 
variety of additional differential equations he may desire to use. We 
note that the functional dependence of |j } on state vector elements is 
not restricted to linear form. In fact, the additional differential 
equations may include sample data systems, hysteresis, dead- 
bands and any other type nonlinearity or discontinuity that can be 
numerically represented. 

The required kinematic and selection transformations, the equations 
of dynamic equilibrium, and the additional controller equations collected 
together are: 

■[■];’((»). 

(*l ■ E[‘]i I'll 

J (12) 

where the subscript j ranges from 1 through the number of cluster bodies^ 
These first order differential equations thus represent the state equa- 
tions governing a system of controlled, interconnected flexible bodies. 
They are clearly of the form 




(13) 


which is easily adaptable to numerical integration. 


Numerical Linearization for Frequency Domain Studies 

If an equilibrium state is known, and if it is desired to examine 
stability characteristics with respect to that state, then numerical 
perturbation techniques can be employed. Considering an autonomous 
system, such that fk = 0; and expanding Equation (13) in a Taylor’s series 
(in Ayk) yields the linear, perturbation equations 


I Ay I = A j |Ay| 


(14) 
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with 


A. . 
ij 


!£i 

ayj 


( 15 ) 


The DISCOS program system obtains the partial derivatives 
numerically (see Reference (10)). Now, as the original system states 
don't necessarily include observed outputs of the "plant", the program 
system also linearizes additional functions that relate observed outputs 
and control actuator torques (forces) to the original states, such that 


' ['=] 


(16) 


The partial derivatives, C^-, having been numerically defined, pro- 
vide a means of eliminating as many of the original perturbated states 
as there are elements of | Aw { , and thus, a replacement similarity 
transformation on the linear equations (Equation (14)) can be executed. 
That is, it is possible to develop 


H '[*]{'} 

such that 

or 

{'} ■ [**]{"}• 


(17) 


(18) 


The linear equations (Equations (18)) represent the exact same 
system as Equation (14), ( [A*] is similar to [a] ), but the new pertur- 
bated states now include observed outputs and actuator torques. The 
advantage of transforming the system's linear equations to the form of 
Equation (18) is that now a plant-controller concept is inherent. Thus, 
subpartitions of [A*] can be easily manipulated so as to create charac- 
teristic matrices corresponding to forward- loop , return- loop and loop- 
gain definitions of the dynamic system. In summary, [A*] contains all 
of the required numerical information to provide the general format of: 

{x} - [aj {x} + [bj {u} (19 

corresponding to any open or closed-loop block diagram of the dynamic 
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system. Thus, there is provision for conducting frequency domain studies 
that support stability analysis and control system design. 


Analysis of a Multibodv Spacecraft 

An analytical investigation utilizing the DISCOS digital simulation 
code was conducted to investigate the d 3 mamical behavior of a typical 
satellite configuration. The configuration chosen was a spinning vehicle 
possessing four doubly hinged solar array assemblies, as shown in Figure 
1. The vehicle was symmetric but consisted of two types of solar array 
subassemblies. Upon deployment, the x-z plane assemblies deployed such 
that “x “ 90° and 7^ - 180°. The y-z assemblies deployed in a similar 
fashion; however, they also articulated such that upon deployment, 

V'y = 90°. 



Figure 1. Spacecraft Schematic 


The analytical investigation was two- fold in that a nominal config- 
uration without the array constraint was first exercised to verify simu- 
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lation of the deployment mechanism. The spacecraft was despun by apply- 
ing a torque directly to the main body and, as the despin occurred, the 
arrays were simultaneously deployed. The array deployment then resulted 
from both spin induced dynamics and the spring/dashpot/cable deployment 
mechanism characteristics. Nondimensionalized results for the nominal 
despin/deployment are presented in Figure 2. The figure also shows a 
comparison with results from the NBOD simulation of Frisch (Reference 4). 
The variation in x and y array deployment characteristics represents the 
time delay of y assembly deployment initiation due to the "yo-yo" despin 
cables uncovering first the x arrays and then the y arrays. 

The constrained array analysis employed a flexible linkage between 
adjacent arrays to affect the desired non- topological tree type config- 
uration. Flexible modal properties were also considered for the panel 
portion of the array assemblies. The spacecraft and the spar portion of 
the arrays were considered as rigid. Results for the despin/deployment 
dynamics of the non- topological configuration are given in Figure 3. The 
for the X assembly arises because of flexibility between the opposing 
X assemblies. There is no restraint between the spacecraft and the x 
assembly. 

A test program was conducted to qualify the analytical results. 
Comparison of test data and analysis was generally very good, however, 
test results indicated more symmetry between x and y assembly positions 
than did analytical results. This was attributed, in part, to tolerance 
accumulation in the structural joints. 


Discussion of Selected Applications 

The system governing equations of motions, as incorporated in the 
DISCOS program system, have been applied to numerical analysis of a wide 
range of dynamic systems. The results thereof have been used as both 
design aids and to confirm ground test and flight data. The following 
text briefly describes selected efforts and is included here to provide 
the analyst with insight as to general applicability. 

Ap pendage Deployment 

The SCATHA (^acecraft Charging At High Altitude) spacecraft will 
deploy five dual element tubular experiment booms in three sequential 
events. The vehicle is spin stabilized; no active control system exists. 

Spacecraft/boom design was accomplished following development of 
the deployment time history; the deployment profile influences stability 
considerations and deployment rates at boom latch were used to size the 
boom tubes and the boom latch mechanism. Additionally, the analyses were 
used to show that the acceleration environment caused by the boom latch 
event would not be detrimental to the experiments. 


99 





^*^ 9/9 uoTljsod 


100 


Figure 2. Nominal Deployment Time Histories 





stagin g 


The SCAIHA mission has the requirement to jettison the Apogee 
Insertion Motor (AIM) following its burnout. The AIM jettison system 
consists of a single-point spring separation device; there are no guide 
rails to constrain motion. The separation analysis of the two- body 
system is further complicated by the fact that the spacecraft is spinning 
with a finite coning angle, neither separating body is inertially symme- 
tric, and the separation device (spring) may have a lateral force com- 
ponent. Bounds on the design coning angle at separation, allowable 
inertial asymmetry and separation device tolerances were established 
to preclude the possibility of short term recontact of the AIM, Addi- 
tionally, an envelope of the separation relative velocities was developed 
for use in long term recontact studies. 

Orbital Docking 

The orbital docking maneuver of the proposed Space Tug was investi- 
gated in order to establish the effects of large amplitude fluid motion 
which might occur during the impact event. A variable length pendulum 
analog was implemented to account for fluid mass motions and the time 
histories of gross vehicle motions and interaction forces were established 
to serve as an aid to impact attenuation mechanism design. 

Digital Control Law 

The Block 5-D spacecraft seirved as the prototype for the TIROS-N 
satellite and was analyzed in order to provide a valid simulation model 
directly applicable to investigation of the TIROS-N on orbit attitude 
dynamics. The configuration consisted of a central body with an articu- 
lated solar array subject to a four channel controller: (1) a three axis 

digital autopilot employing saraple-and-hold and computational delay logic 
and (2) an array drive motor to conq>ensate for spacecraft motions and 
provide a fixed inertial orientation for the solar array. The control 
law (furnished by Goddard Space Flight Center in block diagram format) 
was cast in first order form and interfaced with the spacecraft equations 
of motion. Response to initial attitude errors were determined. 

Other Applications 

The versatility of the DISCOS program system admits its use for a 
wide range of dynamic studies. The nonlinear time domain response analy- 
sis capability has been, or will be used (in addition to those items des- 
cribed previously) to Investigate (1) space fabrication techniques en- 
visioned for future large space structures, (2) launch vehicle- performance 
associated with the Titan III, Stage II post-SECO event, and (3) the 
librational motion of a very large satellite subject to gravity gradient 
excitation and orbital transfer maneuvers. The frequency response analy- 
sis capability has been successfully exercised in the development of 
stability criteria for the Titan III, Stage I POGO event. 
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Concluding Remarks 


A cursory description of nonlinear differential equations, of state 
that govern dynamic behavior of a system of interconnected flexible 
bodies has been given, along with discussion of an available digital 
computer code that synthesizes and numerically solves the equations in 
time and frequency domains. A few of the numerous possible applications 
of the program code are discussed. A comprehensive expose^ of the analy- 
tical and programming considerations is not possible herein, due to space 
limitations and, of course, that has not been the objective of this dis- 
cussion. The objective has been to highlight some of the attendant 
features of the formalism and program code so as to provide insight into 
the motivation and rationale behind them. 

Lagrange multipliers have been used for every constraint of the 
system, not just the ones involved in the closed loops. This has re- 
sulted in increased program flexibility, in that some of the rheonomic 
constraint conditions can be unilateral. 

In conclusion, the authors feel that the nature of the formalism 
and associated program code, DISCOS, is sufficiently general to cover a 
broad range of future applications. 
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Abstrac t 

A unique approach to dynamic response analysis 
of highly flexible, rotating spacecraft is pre- 
sented. Formulation of applicable equations of 
motion includes dynamic coupling between elastic 
.deformations and spacecraft rotation rates. Cou- 
pling is shown to be a closed-loop effect, in which 
elastic deflections significantly influence vehicle 
angular rates, and conversely, elastic response is 
strongly influenced by (if not entirely caused by) 
vehicle spin rates. Angular motion is shown to be 
most strongly influenced by changes in inertial 
characteristics caused by elastic deformation, and 
this deformation is in turn Influenced by centri- 
fugal and Coriolis forces. Equations of motion are 
nonlinear, first-order differential equations of 
first degree. State variables include projections 
of translational and rotational momentum vectors 
onto the rotating vehicle frame and the generalized 
momenta and deflections corresponding to a set of 
normal vibration modes. Additional state variables 
include inertial position coordinates and quater- 
nion (Euler 4-parameter) elements to characterize 
vehicle attitude. Numerical results for a typical 
spacecraft are presented, and the influence of 
flexibility is evaluated. 

Introduction 

A variety of mathematical formulations and com*' 
putatlonal procedures related to dynamic behavior 
of spinning flexible satellites has appeared in re- 
cent literature. Most Investigators have been con- 
cerned with rotational stability of nearly rigid 
vehicles undergoing small elastic deformation. In 
addition to intrinsic dynamical features of a 
freely rotating flexible vehicle, other authors 
have considered the effects of external distur- 
bances such as gravity gradient and attitude con- 
trol . 

Approaches to derivation of governing differ- 
ential equations of motion are diverse. Likins'^^ 
has classified analytical procedures in three cate- 
gories: (1) discrete coordinate methods, (2) hy- 

brid coordinate methods, and (3) vehicle normal co- 
ordinate methods, and has found each to have an 
area of applicability. With his hybrid coordinate 
method, Likins separates a vehicle into a number of 
idealized structural subsystems, each of which is 
classified either as a flexible appendage or as a 
rigid body or particle. Applications are pre- 
(2,3) 

seated to demonstrate the method's utility in 

design of an attitude control sy.stem, and in analy- 
zing a dual-spin spacecraft with solar panels and 
a damped linear oscillator that simulates a nuta- 
tion damper. 
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The Eulerlan N-body dynamic formulation of Hooker- 

Margulies^^^ and Roberson-Wlttenburg^^^ has been 
programmed for general use by Farrell, Newton, and 

Connelly. With this method, the vehicle is con- 
sidered as a set of N interconnected rigid bodies. 

In a recently published two-part paper, Ness and 

Farrenkopf and Ho and Gluck^®^ extend the N-body 
concept with two methods. The first, called the 
unified approach, features an automated synthesis 
of the governing differential equations of motion 
of a multlbodied flexible spacecraft system. With 
this method, terminal bodies of the topological 
tree configuration (no closed loop) may be flexible, 
but interior bodies must be rigid. The second 
method, called the perturbation approach, is more 
general in that any body is allowed to deform, but 
spacecraft motion due to flexibility is treated as 
a small perturbation of nominal motion, resulting 
in piecewise linear differential equations. 

(9) 

Keat * describes a method analogous to the N- 
body approach that is applicable when satellite non- 
rigidity is characterized by generalized coordi- 
nates. Other approaches to derive governing dif- 
ferential equations include a formal application of 

Lagrange's equation by Newton and Farrell . 

Meirovitch^^^^ uses Hamilton's principle, where 
body motion is described by a "hybrid" system of 
equations (not to be confused with the hybrid co- 
ordinates of Likins) consisting of both ordinary 
and partial differential equations. 

This paper presents a new derivation of appli- 
cable differential equations of motion. The ap- 
proach considers a more fundamental basis, namely 
the use of D'Alembert's principle. This approach 

( 12 ) 

is somewhat similar to a method reported by Likins 
in one of his earlier papers. However, the proce- 
dure presented here results in a system of first- 
order ordinary differential equations, akin to 
Hamilton's canonical equations. It also does not 
require use of normal modes of deformation for ei- 
ther the spacecraft as a whole or for appendages as 
separate structural components, but may use normal 
modes to facilitate numerical solutions. 

The method has the following advantages: (1) 

the final format of governing differential equations 
is similar to classical Newton-Euler rigid-body 
equations, (2) derivation of equations is uncompli- 
cated, (3) vehicle strain energy, kinetic energy, 
and total angular momentum can be easily monitored, 

(4) equations are particularly amenable to program- 
ming for digital simulation. 

We have limited the scope of this study to ex- 
clude effects of external environment, being con- 
tent to examine the behavior of a freely spinning 
flexible satellite in the time domain. Extension 
of the analysis to consider attitude control, gra- 
vitational forces, and other disturbances presents 
no great difficulty. 
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Development of the Analysis 

Governing nonlinear differential equations of 
motion for spinning flexible spacecraft result 
from direct application of D'Alembert’s principle. 

( 13'i 

Lanczos'- ' states that all the different prin- 
ciples of mechanics are merely mathematically dif- 
ferent formulations of D'Alembert's principle." 

He implies that D'Alembert's principle is a truly 
fundamental principle, particularly useful when 
using "kinematical variables" co characterize the 
motion. 

In studies of rigid-body dynamics, it is con- 
siderably more convenient to use angular velocity 
components that are projections of the spin vector 
onto the moving axes. We believe it is more con- 
venient to use kinematical variables when Cvonsider- 
ing the dynamics of a flexible spinning spacecraft. 

Geometry and Kinematics 

Consider the geometry of a deformable spinning 
spacecraft, as shown in Figure 1. 



j Frame 

Figure 1 Geometry of Deformable Spinning 
Spacecraft 

We use two rectangular Cartesian frames of 
reference; an inertially fixed (Newtonian) axis 
system (X,Y,Z) and a moving frame (x,y,z). We 
should emphasize that reference point R does not 
necessarily coincide with the mass center of the 
deformable body, nor is it necessarily fixed to 
some material point of the body, although it could 
be either. We represent body deformation as a 
single-valued displacement field n to be measured 
in the moving reference frame. In the figure, 
material point p is shown in a deformed configura- 
tion^^specif ied by the instantaneous position vec- 
tor p. The corresponding point p' is in the un- 
deformed configuration and is specified by the po- 
sition vector p fixed to the moving reference 
frame. These geometric notions are in line with 
traditional treatment of a system of parti- 
cles, although our development considers a 
continuous structure. 

In view of these geometric considerations, we 
may now specify an absolute velocity field, v, 
associated with the deformable body, and by kine- 
matic principles, express velocity of the typical 
point p as 

v-v^ + ;:;xp + (p>^ 


where v denotes the absolute velocity of refer- 
R 

ence point R, w represents angular velocity of the 
moving reference frame, and (*)j. represents the 

time derivative of the enclosed quantity, as seen 
by an observer rotating with the moving reference 
frame. In particular, in Eq 1 the term (p)^ can 

be written as 

(p)f = (n)j. ° ^ \ ^ '^y ^ 

The displacement field, n, can be further de- 
veloped as a linear combination of a finite number 
of single-valued displacement functions (of x, y, 
and z) , so that we can write 

n (x,y,z,t) = 

k 

and it follows that 

(n)^ = £ \ (n.y.z) . (4) 

k 

We can also express the acceleration field, a, 
associated with the deformable body as 

a = (v)^ + Qj X V (5) 

which will be needed in subsequent expressions of 
equilibrium. We have tacitly implied that compo- 
nents of all vector quantities (except r^^, the po- 
sition vector to reference point R of Fig. 1) are 
projections onto the rotating reference frame (i, 
j, k) , thus we can write 

v^ = iu + jv + kw I 

ti) = io) + ju) + ko) I 

x y z I 

p = I (x + E «k) + ^ (y + S V 0 I 

(z + E^k ^k) • ) 

With the preceding definitions and notations, we 
are in a position to develop applicable state 
equations . 

Development of State Equations by Direct Applica- 
tion of D'Alembert's Principle 

Equations of motion for any dynamical system 
are derived, in one way or another, from consider- 
ations of equilibrium; some authors use the term 
equations of equilibrium synonymously with equa- 
tions of motion. A system of impressed forces is 
not generally in equilibrium, and motion of the 
system (implying accelerations and additional in- 
ertial forces) results in an overall system of 
forces (impressed and inertial) that is in balance. 
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D' Alembert's principle states that any system 
of impressed and inertial forces is in equilibrium, 
and therefore, the total virtual work of the ays* 
tern is zero for reversible virtual displacements 
that satisfy the given klnematical conditions, or 


6r • (f - oa) dV 


D'Alembert's principle, restated as Eq 7, is the 
basis from which we derive applicable equations 
of state. 

The given kinematical conditions with which 
the virtual displacements must be consistent are 
represented by the velocity expression of Eq 1. 
Elements of this expression are kinematical vari- 
ables (sometimes referred to as nonholondmlc ve- 
locities) because we stipulated chat all vector 
quantities would be represented as projections on- 
to the rotating axis system. These nonholonomic 

velocities (u, v, w, w ", u , m ) are not time de- 
X y 2 

rivatives of actual position coordinates. Multi- 
plying the nonholonomic velocities by differencial 
time (dt) leaves a sec of Infinitesimal quantities 

(dr , dr , dr , d0 , d0 , d0 ) that are not dif- 
X y 2 X y 2 

ferentlals, but are quantities referred to as 

"differencials of quasi coordinates". 

Because they are infinitesimal quantities, they 
correspond to a sec of virtual infinitesimal dis- 
placements (6r , 5r , 6c , 60 , 66 , 50 ) and, in 
X y 2 X y 2 

view of Eq 1, 6 and 7, we can write the virtual 
displacement field as 


= l( 6 r +60 p - 66 p + ^ . 6 C 1 

\ X y z 2 y k y 

+ j(6r + 50 p - 50 p + ^ ^ 6^,1 

\ y 2 X X 2 ^ 

+ k(6r +60 p - 60 p + ^ y it 6^, ) 

\ z X y y X ^ zV. k/ 


Combining Eq 7 and 8 results in 

f |(v) + ;; X V I odv = /" f dv 

•'v •'v 

J P ^ (v) + oj X v odV = / D X f 

»» * r * •'it 


p ■ / 0 * ’ 

“ -'v 

’’j ' f \ 

^Ic V 


0 * V odV 


• v odV 


If we differentiate left and right members of 
Eq 12 through 14 and compare the results with 
Eq 9 through 11, we arrive at the following state 
equations 


(0 ■ X 


pxfdV-tDXp -v„xp 
w R V 




f dV - U» X p 


/ (K 


We see immediately chat, if are suppressed, 

the first two of Eq 15 represent the rigid-body 
equations of motion. The third of Eq 15 is, ne- 
glecting Che second term on the right, a familiar 
equation of motion for the kth normal coordinate 

if we consider ' f dV to include effects of 

structural damping and restoring forces and if ve 
define to be a normal coordinate. 

Referring to Eq 12 through 15, and using the 
momentum components defined by Eq 16, viz., 

p.. “ h\ + jy2 + ^y3 


p„ = iyu + Jys + kyb 




dV ( 10) 


allows us to write the state equations as 

Vl = Gi - (a> ,3 - - (vy^ - wy^) 


^2 “ S • ■ V 3 ^ ■ ^"^4 - “V 


(v) + li) X V odV 


because- virtual displacements 6r , 6r 6£. 

X y ’ k 

are arbitrary and independent quantities. These 
results are familiar expressions of equilibrium 
if we recall some results from studies of dynamics 
of a System of particles. 

Let us now define Che following momenta 


yj = Gj - (ui^y^ - - (uyj - vy^) 

^4 = ^=4 - <“y^6 - “.^5^ 

yj = G 3 - (.^y^ - 

^6 ° • V 4 ^ 

^k +6 ° S +5 • 


It can now be seen that the (6+n) state equa- 
tion given by Eq 17 are expressed as functions of 
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the state variables themselves and other undeter- 
mined quantities. It Is therefore necessary to 
develop a means to establish these quantities in 
terms of known state variables. 


If we substitute Eq 1 In Eq 12 through 14, 
respectively, there results 


Pv = I (- 


+ (1} X p + 




odV 


P X Vr + p X (oi X p) 


V 

+ * k') « 


'•'k' ‘■k 


odV 




■l\ 


♦k ■ ■ (P ^ 


+ E ♦k • s 


OdV 


(18) 


(19) 


( 20 ) 


Integrating the first of these (over the volume V) 
term-by-term allows us to write 


p “ m V. + u « 
•^v R 


^?odV.E(//, odv)j, 


“TnVR + tii><S + 



(21) 


where m is the total system mass, S is the 
"static" moment vector. 

In a similar manner, the other equations 
yield 

Pw = ® h +Edk (22) 

k 

Pc^.= \ • \ + ^ 

where h is a momentum vector referred to the mov- 
ing coordinate system and is defined as 


h * li j kj (J) {(*)} 


(24) 


and d, and 
k 

Comparison 
the matrix 
ables (y^^, 

ables (u , 

X 


e. are undetermined coefficients. 

kj 

of Eq 21 through 23 with Eq 16 yields 
relationship between the state vari- 

yo****» y^Li ) the kinematic vari- 
2 b+X 


' ’^6+k^ 


’^l 


J -J -J 

XX xy xz 


-s 

z 


d , . 

xl 


. , d " 
xn 

^2 


J -J 

yy yz 

S 

z 


-S 

X 

d , . 

yi 


d 

yn 

^3 


J 

zz 

-S 

y 

S 

x 


d , . 

zl 


d 

zn 




m 



a , • • 

xl 


. a 

xn 

5'5 

= 



m 


"yi • 


. a 

yn 

^6 





IR 

Si • 


. a 

zn 



(Syrametri 

;) 



^1 * 
®22 


• ®i 

In 

/6+n_ 








e 

nn 




(25) 


where the coefficients in the above array are 
defined to be 


J = (J 
aB ^ aB 


>o" Efosj" Vj) 


J 

aa 


(J ) 
aa c 


+ 2 



+ b 




^ S '^ajBk S ^k 


n n ^ 

+ E E ' 

J k \ 


“^Sjek * Sjyk)S 


/ a = X 
when ! (1 * y 


5 « y 

6 = 2 
6 = X 


Y = z 

Y = X 

Y = y 


(26) where a, 6 ® x,y,z and a 6 

\ - E Ej, 

where a = x,y,z ; 


(27) 


( 28 ) 
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^ak “ 

’’Brk ■ 

'vBk * 

n t 

^ (‘^BjYk ■ ‘"yjBk) 

(29) 


/ O « X 

6 = y 

Y e 2 


when 

' a = y 

B = z 

Y = X ; 



^ a = z 

6 = X 

Y = y 


"jk = 

"^xjxk ^ 

“^yjyk 

+ c 

zjzk 

(30) 

where 

j.k = 1. 

2,--- 

. o ; 



and 



(31) 


In the definitions above, we have used 

(J „) = / (6^ + Y^) adV 

““ ° -'v 


1“ ° 

X 

B = y 

Y ' 

z 

when ! a » 

y 

6 = z 

Y = 

X 

(a = 

z 

B = X 

Y = 

y 


(J „) = / aSodV 

ae o 


where o,6 = x,y,z and a i* 6 ; 



where a = x,y, r ; 



(32) 


(33) 


(3A) 


(35) 


Remaining n forces acting In deformation coordin- 
ates, C, are defined to be 


^k+6 “ ®k+6 \+5 n 


(37) 


where 




®k+6 = \ • f dV 


(38) 


and 


b , + b , + ^(c , , + c , .) C- 

yyk 2 zk V ykyj zkzj 


+ “y '’zzk + ”xxk ^ ^^“=zkzj + <=xUxj> 


'’xxk * **yyk * ? ^“^xkxj ^ *^ykyj^ 


- u w 
X y 


o CO 0) 1 

X z 


^xyk ^ ^yxk ^ ^ ^*^xkyj ^ ^ykxj^ 
^xzk ^ ^zxk ^ ^ ^'^xkzj ^ ^zkxj^ ^j] 


- (0 to 

y 2 


b ,+b ,+5irf(c, 

yzk zyk ^ ykzj zkyj j 


+ to [wa , - va , 1 + to [ua , - wa , ] 
x yk zk' y zk xk" 


+ mjva^^ - 


“x ^ (Cy^zj “^zkyj^ 


■*• “v S (c,,.„^ - c„^,,) 


zkxj xkzj j 


“z ^ ^"^xkyj ■ 'ykxj^ ■ 


(39) 


where a,S * x,y,z; j,k = n . 

Solution for unknown velocities on the right-hand 
side of Eq 25 in terms of state variables is now 
straightforward. 


The first six forces (torques) in the state 
equations are components of the resultant im- 
pressed force (torque), namely 


G = y (p f) dV 
V 

‘=1+3 


1 “ 1.2,3 


(36) 


All previously defined coefficients depend on 
the variables, and therefore, the state equa- 
tions must be augmented to establish values for 
the respective This can be done by integra- 

ting the which are obtained from solution of 

the simultaneous equations given by £q 23. The 
final form of the (2n+6) state equations Is 

Yj = fj (y^. yj.---. y2„+6- O • (W 
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This presentation of the state equations is 
complete in that all variables required to char- 
acterize vehicle elastic response have been in- 
cluded. However, it is often necessary, or at 
least desirable, to monitor the vehicle's posi- 
tion and attitude with respect to the inertially 
fixed reference frame. This can be done by ad- 
dition of supplementary state variables. The 
position of reference point R can be monitored 
through definition of the three projections of the 
vector onto the inertial axes. These can be 

obtained by integrating components of v^^ along 

the inertial axes. The inertial attitude of the 
moving reference frame can be monitored by using 
quaternion (Euler 4-parameter) elements. We have 
found this to be more practical than use of Euler • 
angles because the potential of "girable lock" as- 
sociated with large excursions is eliminated. 

Application 

Equations of motion previously described were 
programmed for solution on a digital computer. 

The program was then used to analyze the dynamic 
response of a representative elastic spacecraft, 
the Phase A version of the proposed Planetary 
Explorer (Fig. 2). This investigation was in- 
tended to establish the influence of variations 
in flexibility on overall response. 



Figure 2 Planetary Explorer Venus Multi-Probe 
Configuration 

Dynamic Model and Response Analyses 

Three response simulations were performed. 

The first, Case 1, considered the body to be 
rigid, and the other two considered variations in 
overall structural flexibility. Case 2 assumed 
nominal flexibility derived from a finite-element 
representation of the structure, and Case 3 assumed 
a 25% increase in flexibility. Modal character- 
istics required as input to the response analyses 
were obtained from the 60-degree-of-freedom model 
shown in Figure 3. The first six elastic modes 
are shown in Figure 4. 




The structure was assumed to be free from the 
influence of external forces (torques); elastic 
response was excited from the undeformed state 
through Imposition of initial reference-axis angu- 
lar velocity components 


w » ll j kj 


0 

I 

■n , 
lOn 


(rad/sec) 


Because normal vibration modes were used, and 
reference point R was initially at the center of 

mass, vectors a, and S are zero (a, - 0 from 
k o k 
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orthogonalicy considerations of the unsupported 
vibration modes). Therefore, the mass center and 
reference point remain coincident throughout the 
simulation. 

Results 


The three principal inertias and projections 
of the spin vector onto the principal axes are 
shown in Figures 5a through 5c and Figures 5d 
through 5f, respectively. These inertias, as a 
function of time, are determined by extraction of 
characteristic values (eigenvalues) of the inertia 
dyadic, J, of Eq 2h and 25. Here we note with no 
great amazement that, for the rigid case, the 
principal inertias are constant, and spin vector 
components on the intermediate and minor princi- 
pal axes are oscillatory; one trace leading the 
other by 90®. The frequency of these oscillations 
agrees with that predicted from the rigid-body 

. ( 14 ) 

expression 




V 




- J ) (J - J ) 
XX 22 yy 


J J 
XX yy 


obtained by considering the stability of pertur- 
bations in 0 ). 

Examination of the principal Inertias for the 
two flexible cases again indicates an oscillatory 
nature, the steady-state value of the major axis 
Inertia being somewhat increased over the initial 
value. This is consistent with the characteris- 
tics of the spin vector components, where we see 
that the steady-state projeccion on the major 
axis is decreased from the initial value. This 
is a direct result of conservation of angular mo- 
mentum. We note further chat the projections on 
the intermediate and minor axes indicate a null 
value at steady state. It can therefore be con- 
cluded that, in the presence of flexibility (and 
implied structural damping), the spin vector tends 
to align with the major principal inertia axis. 

The orientation time history of the z- 
reference axis (k of Fig. 1) to the principal 
axis system is shown (Fig. 6) for the two flexi- 
ble cases. These plots indicate angles determined 
from three direction cosines extracted from the 
orthogonal rotation transformation that relates 
reference and principal axis systems. 

For nominal flexibility (Case 2) , we note 
that the reference axis is nearly coincident with 
the principal axis system and that steady-state 
excursion will be less than 1®. Case 3 shows 
markedly larger excursions, for here we note that 
although one angle is approximately 90®, angles 
between the reference axis and the other two prin- 
cipal axes have a significant value. In fact, the 
steady-state excursion is seen to be about 10®. 
This increase in reference-axis excursion results 
from increased structural flexibility, although it 
is not possible to specify the exact relationship 
between flexibility and reference-axis excursion. 

Figure 7 shows three characteristic motion 
variables for the two flexible cases. Figures 7a 
and 7b indicate the projection of the spin vector 
onto the momentum vector, which, as we prescribed 


no external torques, is known to be fixed in in- 
ertial space. This projection is oscillatory and 
approaches a steady-state value whose magnitude 
is somewhat less than the initial value, implying 
a steady spin condition. Inspection of the 
steady-state value and comparison with Figures 5e 
and 5f, where the spin vector tends toward align- 
ment with the axis of maximum principal inertia, 
allows us to conclude that the major principal 
axis tends to coalesce with the momentum vector. 
This conclusion is verified through Inspection of 
Figures 7c and 7d, where we have shown the time 
history of the angle between the spin and momentum 
vectors. Note chat the angle is a measure of the 
nutation and that structural flexibility (damping) 
acts as a nutation damper; the damping Increasing 
with increasing activity in the deformation 
coordinates 

The angle between the z-reference axis and 
2-lnercial axis (K of Fig. 1) is shown in Figures 
7e and 7f. These data indicate the relative ex- 
cursion of the spacecraft with respect to the dn- 
ertial frame. 

The elastic displacements at the three satel- 
lite points (the probes in Fig. 2 and 3) are 
shown in Figure 8. For nominal (Case 2) flexi- 
bility, the magnitudes of the three n deflections 

z 

are approximately equal and the n and n deflec- 

X y 

tions indicate the almost symmetrical axial exten- 
sion of the booms. This result is not surprising, 
for the Imposed initial conditions resulted in a 
spin rate that was in close agreement with ob- 
served natural frequency of the third elastic mode, 
and the effect of this spin race has been to ex- 
cite an almost pure mode. No such observation is 
evident for the modified flexibility case because 
the resultant system natural frequencies are much 
lower in value than the imposed spin rate. In- 
creased structural deformation with increased 
flexibility is, however, quite evident. 

Figure 9 shows time variations of the flexible 
system energies, where we noted the constant 
rigid-body value on the kinetic-energy history. 
Kinetic energy (Fig. 9a & 9b) is seen to decrease 
and approach a steady-state value. This further 
validates previous assertions with respect to 
variations in principal inertias and angular ve- 
locity components. The steady-state value can be 
verified through inspection of Figure 5. Figures 
9c and 9d indicate the corresponding strain en- 
ergy arising from structural deformation. Once 
again we see the tendency toward the steady-state 
value that will be achieved as elastic motion is 
finally dissipated. Case 3 indicates the higher 
value of strain energy expected because of the 
nature of elastic deflection histories shown in 
Figure 8. 

In conclusion, we point out that total system 
energy (Fig. 9e & 9f) is monotonica.l ly decreasing 
to a value equal to the sum of the steady-state 
kinetic and strain energies. 

The results — that the spin vector coalesces 
with the axis of maximum principal inertia — are 
in agreement with the precept that a spinning sys- 
tem will, in the presence of damping, seek the 
lowest possible energy level consistent with a 
constant angular momentum. 
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Conclusions 

We have used a direct application of 
D'Alembert's principle to arrive at canonical 
state equations dnd have shown that these equa- 
tions have minimal coupling. The similarity of 
this development to classical Newton-Euler equa- 
tions for a rigid body clarifies assessment of 
effects of deformation on the overall motion. 

This clarity is not so evident when the equations 
are derived by application of Lagrange's equations 
or Hamilton's principle. 

Numerical examples yield results that agree 
with previous knowledge regarding behavior of 
flexible spinning vehicles. In particular, 
coalesence (in the steady-state configuration) of 
the spin vector, momentum vector, and major prin- 
cipal inertia axis was observed. Total system 
energy was shown to decay to a steady-state value 
consistent with a constant angular momentum. Re- 
sults clearly show that Initial and final princi- 
pal axes can be misaligned depending on relative 
values of spin rate and structural vibration fre- 
quencies. The steady-state configuration is not 
necessarily axisymmetric, even chough the unde- 
formed configuration was. 

We have purposely limited the numerical in- 
vestigation by excluding effects of gravity or any 
other external disturbances. In so doing, we 
have been able to isolate mutual effects of vehi- 
cle spin and deformation; in themselves complex 
phenomena. 

By concentrating on a freely rotating flexible 
vehicle we have established a practical and mathe- 
matically rational procedure for analyzing more 
complex dynamical systems. In particular, exten- 
sion of the methodology to include gravity gradi- 
ent, attitude control, or dynamics of an inter- 
connected system of flexible vehicles (e.g., ar- 
ticulating appendages, multispin configurations, 
etc) can be accomplished in an orderly fashion. 
Another extension might be linearization of the 
equations to examine stability in the frequency 
domain. 

In conclusion, we note that the procedure can 
consider a rigid main body and, therefore, might 
be used to support a test program in which one or 
more flexible structures is attached to a rigid 
spin platform. 
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APPENDIX A 


DEVELOPMENT OF THE INERTIAL INTEGRALS 


In developing the equations of motion (paragraphs II.B and II. D of the main text), certain 
inertial integrals are identified that are required to account for the deformation-dependent 
inertia matrix and that are involved in calculating the effects of centrifugal and Coriolis 
forces. 

The basis for calculating these integrals is a triple-matrix product that involves a so-called 
discrete mass matrix, [M] , which is assembled by finite element techniques and which may 
be used in calculating vibration modes. The other constituent of the triple-matrix product 
is a modal transformation that transforms ordinary velocities associated with the finite- 
element model to the velocities of the {u}^ vector. 

Referring to the transformation as [0] , the triple-matrix product is therefore 


[m] = [ 0 ]T [M] (^] (A-1) 

which is the basis of the kinetic-energy expression of equation 11-21 of the main text. Now, 
the mass matrix, [M] , is invariable with respect to the body’s deformation. However, the 
modal transformation, [ 0 ] , depends on the {^} in a linear fashion, or [ 0 ] may be expanded 
as 


[<^] = [</’lo (^'2) 

where [0]^^ is a matrix of constant elements, and [A0] is variable with respect to deformation. 

On substituting equation A-2 into equation A-1 and referring to equation 11-86 of the main 
text, it follows that 


K1 = [M] [0^] 

(A-3) 

= [A0]T [M] [0J + [0JT [M] [A0] 

(A-4) 


Preceding P^ge Blank 
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and 


KljkMk = M [A«] (A-5) 

Assume that the finite-element model of the body has a “global” Cartesian frame in which 
the ordinary velocities are measured, and further assume that the generalized coordinates 
of the finite-elementi model are grouped (or ordered) so that all the x-translations are together, 
followed by all the y- and then z-translations and that the translations are followed by sets 
of X-, y-, and z-rotations. With this implied ordering, it follows that the discrete mass matrix 
is partitioned in the form; 


(A-6) 


with p, q, and r corresponding to rotation coordinates about x-, y-, and z-axes, respectively. 
Similarly, the iriodal transformation is partitioned as 


[M] = 


m 


m 


xy 


m 


m 


xp 


m 


xq 


m 


m.. 


yy 


m. 


yz 


m.. 


yp 


m 


yq 


m.. 


yr 




tn 


zp 




zq 


m 




pp 


m 


pq 


m_ 


pr 


(symmetric) 




qq 


m. 


qr 


m. 


[ 0 ] 





{If 




i^4 


{’1+14 


M 


[h,l 

|y^y} 





M 

[hj 

{if 






KJ 


Ilf 








{If 






(A-7) 


Each square subpartition of equation A-6 has rows equal to the number of structural joints 
(collocation points) of the finite-element model, as does each subpartition of equation A-7. 


A-4 







The submatrices in the last column partition of equation A-7, ([h^ ] , [h^ ] , have 

columns equal to the number of deformation modes used to represent the body and are 
matrices of modal translation and rotation amplitudes. 

The forms of [0^ ] and [A0] are seen immediately from equation A-7 in that the only 
nonzero parts of [A0] are attributable to the {i?} vectors. The [0] matrix is effectively 
a kinematic velocity transformation consistent with the form of equation IT25 of the main 
text, and it follows that 


{’’xl = (M (Sl 

(l,t * |1>,1 jSl (A-8) 

K( = W |fl 


Equation A-4 shows the product of two constant matrices; namely, [M] [0^ ] . The two 
triple products on the right of equation A-4 require evaluation of only the first three row 
partitions of [M] [0 q ] . Thus, the following is defined 

(The first three row partitions of [M] [0^^] ) = 



{>’,4 l'’„l 


(Pxal 


Kkl 

I'’,.} 


1%( 



Pykl 

{'’zlf 

{‘■.at 


tf.sl 


K1 


with, for example. 


(A-9) 


= Kp] {U Kzl {y} - Kyi M 


and 


Pxkl = Kxl Kl + Kyi Kl Kzl Ihj 
+ Kp] K1 Kql K1 + [™xrl [®zl 


(A-ll) 


It is unnecessary to expand each partition of equation A-9; the partial product is numerically 
obtained, and the examples of equations A-10 and A-ll are for purposes of illustration only. 
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Now, with reference to the intennediate constant matrices given in equation A-9 and the 
definitions of equations A-4 and A-5, the following inertial integrals* are developed: 


K1 = ' PysJ fhj (A-13) 

{«3l = Kj Pv! - l\3 (A-14) 

K1 = (A-IS) 

tsl = - P.sl (A-16) 

[^1 = {fj [h,] - {P,,l [hj (A-17) 

la,l = (Py4l ' tPx4l fhyl <A-18) 

(«g] = IPysl IK^ - Pxsl Ihyl 

[%] = [Pyel Ih,] - IP^,] Ihy] (A-20) 

[by] = [P,J Ihy] - [Py,3 [hj (A-21) 

[b,l = [P,,3 Ih,} - [P,,3 IhJ (A-22) 

[bjl = [Py3l [\] - {P,3l Ihyl (A*23) 


♦The reader is advised to refer to paragraph H.D, particularly equations I1'88 and II-89,,of the main text. 
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[bj = P,,] [hj - [\] + P,,] [hj - [hy] (A-24) 

[bs] = PxJ - Pyil Pxl + Pya] P.l - P. 3 ] Py] (A-25) 

Pel = Pxal Py] - Pyal Pxl + Pza] Pxl ' Pxa^ Pzl (A-26) 

[CyJ = [h,]^ P,J - [h^ PyJ (A-27) 

[C«] = P.1" PxkI - Pxl^ Pzk] (A-28) 

[C,,] = [hj^ P,J - [hy]^ P,J (A-29) 

[C„] = [h [mj [h ] - [hy]^ [m^] [hj 

(A-30) 

- PJ^ KJ [hy] + \hy [m^y] [hj 


= Pzl" Pzl - Pzl" Kzl Pxl 

(A-31) 

- [h,]^ [m^] [hj + [hJ^KJ [hj 
Paal = Pxl^ l"»yyl Pxl - Pxl^ [-"yxl Pyl 

(A-32) 

- [h,]^ K,] [hj + [h^]^ K,] [h^] 

[C,^] = - P,]^ ImJ [hj + [h [mj [hj 

(A-33) 

+ PJ^ [m^] P.] - P.]^ K.] PJ 


A-7 



lC,3l = 


(A-34) 


+ Ih,]^ [m^l [hj - [hj" [hy] 

- [hj" KyJ f^xl ^ Kx3 
+ [hj^ [m^] [hj - [hj^ [m^] [h^] 


= 


(A-35) 
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APPENDIX B 


DEVELOPMENT OF ROTATION TRANSFORMATIONS 


In terms of Euler angles, the analyst may choose from twelve possible orthonormal rotation 
transformations in order to orient one orthogonal triad with respect to another. For each 
of these orthonormal rotation transformations, there is an associated rotation transformation 
that is not orthonormal and that is used to transform angular velocity projections (onto a 
nonorthogonal vector basis), which are time derivatives of Euler angles, to projections (onto 
an orthogonal vector basis) that are commonly referred to as time derivatives of angular 
quasi-coordinates (cj^, cjy, and co^). 

For digital computation, these transformations can be generated automatically, given a 
selected order of rotation. The purpose of this appendix is to indicate the steps and numerical 
manipulations that are required. To this end, consider one of the twelve types (e.g., a 2-3-2 
permutation) as an illustrative example. 

Consider the two orthogonal vector bases, whose relative orientation we want to describe, 
to be 


and 




(B-1) 


(B-2) 


Now, if , ^ 2 ’ three successive Euler rotations about axes (2-3-2), respectively, 

it then follows that 

W-lT.HVj (B-3) 
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(B-3) continued 


and 


COS0J 


0 

-sinfl j 


0 sindj 


r' 

1 0 


-rf 

J 

0 COS0J 


k' 





i«1 • [T,l {el 


- 


“ 

COS03 0 sinfij 


i 

0 1 0 


T 

0 COS63 


k 





On combining equations B-3, B-4, and B-5, there results 

)a( = [TJ [TJ [T3] |I} 


(B-4) 


(B-5) 


(B- 6 ) 


Now, a 2-3-2 permutation means that the first rotation, (6^), is about the second axis of the 
{ a } basis, the second rotation, (d^), is about the third axis of the { e '} basis, and the third 
rotation, (^ 3 ), is about the second axis of the {e "} basis. 

Consider the correlation between Euler rotations and the corresponding axis (table B-1). 

It is now clear that the elementary rotation transformations, ([Tj ], [T 2 ], and [Tj ]), 
always involve 0 ^ , 6^, and dj , respectively, but any of them may have three different forms, 
depending on the axis associated with its rotation. That is, when 0j (i = 1 , 2, 3) is about 



Table B-1 

Correlation of Euler Rotations and Axes 


Type 

1 

2 

0, about 

1,1 

1,1 

02 about 

SilSI 

0j about 

3,k", 

1,1" 




10 

11 

12 

3,K 

3,K 

3,K 

i,r 

2,1' 


3JT" 

m 

i,r" 


axis (1), then 


[T,] = ri 0 


0 COS0. -sinflj 


0 sinflj COS0J 


when flj is about axis (2), 


[Tj] = COS0J 0 sin0j 
0 ' 1 0 


-sin0j 0 cos0j 


and, finally, when is about axis (3), 


COS0. 

-sin0 

1 


sin0j 

COS0J 

0 

0 


Thus, it is evident that, to create the required orthonormal rotation transformation, equation 
B-6, only a rotation type (table B-1) and the three Euler rotations need to be specified. 
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The associated rotational velocity transformations are developed as follows. Consider, again, 
the 2-3-2 permutation. For this case, it is possible to express the angular velocity vector c3 
in two ways; 


CO 



+ kco 

Z 


and 



On combining equations B-10 and B-1 1, there results 

|co| = [tt] 


or 


• 


CO 

X 


CO 

y 

II 

CO 

z 



COS0, 

0 

sin0. 


0 -sin0j 

1 0 

0 co&6^ 


sin0„ 


cosO„ 




w •• 

0 

0 



0 1 


• 

1 jo 


• 


(B-10) 


(B-11) 


(B-1 2) 


or 


H = 1X3]^ [A]-^ 


(B-1 3) 


Now, the inverse transformation of equation B-1 2 is required for hinge kinematics applications, 
or it is necessary to express 


N ‘ ' [13] 


^[E]-« [E] [A]^y' [T3] 
^[EI [Aji-y^ [E] [T3] 


(B-14) 
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with [E] , an elementary row interchange transformation that, for the (2-3-2) example, is 


[E] = 


1 0 0 
0 0 1 
0 1 0 


and causes 


([E] [A]T) 


to be of the form: 


[E] [A]'^ 


a 0 0 
0 1 0 
P 0 1 


so that 



1 /a 0 0 

0 1 0 
-Pfa , 0 1 


(B-15) 


(B-16) 


(B-17) 


with a = sin 02 and P = cosflj • 

The purpose of introducing [E] was to make the form of equation B-1 7 the same for all 
twelve types of Euler rotations. This is convenient with respect to programming considera- 
tions. It follows that: 


For types 1, 5, and 9, a = cos 02 >| 3 =sin 02 - 

For types 2, 6 , and lOj a=sin 02 . /J=cos 02 - 

For types 3, 7, and 11, a = -sin 02 . P - cos 02 • 

For types 4, 8 , and 12 , o = cos02> P = sin02- 


Also, for each of the twelve types, there is an elementary row interchange transformation, 
[E] , that can be constructed from simple inspection of the permutation integers of table 
B-1 (e.g., 2-3-2). In fact, it is not necessary to actually construct [E] because information 
to construct it is merely applied to [Tj ] (interchanging its rows), which produces [E] [T 3 ] . 
Thus, the velocity transformation of equation B-1 4 can be created for any of the twelve 
possible types with comparative ease. 


B-7 
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TIME DERIVATIVES OF KINEMATIC COEFFICIENTS 


The 1 formulation and numerical implementation of motion equations for the system of inter- 
connected bodies involves a vector of Lagrange multipliers, •{ X|. (Refer to equations ll-l 
and II-6 of the main text.) In order to numerically evaluate | X[, time derivatives of kinematic 
coefficients (velocity transformations) associated with hinges must be calculated. 

With reference to paragraph II.C of the main text, note that, for each hinge, there is a [bp ] 
and a [bq ] matrix of kinematic coefficients. The basic form of these matrices is repeated 
here, followed by the sequence of steps necessary for developing their time derivatives. 

The [bp] array is 


and 


[bp] = - 


[^] 


-1 


[ R ] 

m-" 


[ 0 ] 


[q»m] [%] 


[ R ] [S(‘">]Ti R H [ R ] [h ] 

*-p *^p m** *-p m-**’p'* 

L I I 


[b,] 


[,R„] 1 [0] [RJ [aj 

yj [Sffl I [pR„] I IpR„] [hj 


• • 

To develop [bp ] and [b^ ] , it is necessarj^ to expand the following as 


(C-l) 


(C-2) 


£ 

dt 


( 


[n] 


-1 


[ R ] 


) 


[n-^] 


[ R ] 

‘•q m J 




([q« 


R ] [ R ] 

p-* *^p m-* 


[ R 1 [ 

iq pJ ip 


Rl) 

m ■* I 


(C-3) 
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(C4) 


(l'l ‘ ■ l*'l l,R„l * W [,R„I 

(C-5) 

and 

t PS’l) ■ [p«J " IpM PS’l 

(C-6) 

The 3 by 3 matrix time derivatives defined by equations C-3 through C-6 have factors (also 
3 by 3 matrix time derivatives) that are expanded in terms of previously defined quantities 
as follows; 


!pR„l - PK* (lpR„l (%i \U) 1 IpRp,) 

(C-7) 


(,RJ - 1SK-([,R„1 1P,1 IIJ) )[,R„1 

(C-8) 


ip*,i = ipR.i 

(C-9) 

with 

[S2<P)1 -SK*^[^R„| (jppj + Ppl (L!) 

- lpR,l W (Kl * KHU)) 

(C-10) 


IpRpl - lpR,l l,R.l + lpR,l [,R„1 

(C-11) 


(SL”p’l = [sK' ( [l>pl is„l)] 

(C-12) 
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(C-13) 


= [sK* ( lh,l J)] 


Finally, the time derivative of requires additional consideration. (Refer to Appendix 
B.) The rotation transformation, , is developed as 


= 



(C-14) 


and it is shown that the form 


([E] [A] 



[A] 


1/a 0 0 

0 1 0 

-^/a.O 1 


(C-15) 


holds for each of the twelve possible Euler rotations. In that [E] is constant, [A] depends 
only on 62, and [Tj ] depends on 6^ only. It follows that 

w = ^ w m ff,) + K (w [El) [T,i <c-i6) 

• • 

where the Euler angle rates, 62 and , are numerically evaluated before their use in equation 
C-16 through application of equation II-3 of the main text (that is, they reside in that part 
of the state-vector time derivative, | y| , that has been evaluated). 
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MONITOR OF SYSTEM MOMENTA AND ENERGIES 


Developing state equations for predicting dynamic response of a system of interconnected 
flexible bodies involves a considerable amount of compUcated formulation and programming 
codes regardless of the particular method of analytical mechanics selected for basing develop- 
ment. 

The inherent complexity of such a digital simulation program gives rise to the question: Is 
there any way to check the program validity? In an attempt to answer this question, one 
suggestion is to compare results with those of other dynamic simulations or hardware tests. 

If such a comparison is positive, credibility (to a degree) is established. However, another 
absolutely necessary (if not sufficient) condition must be passed to establish validity. For 
a dynamic system free of external forces and torques, angular and linear momenta must be 
conserved, and total energy (kinetic plus potential) must not increase in time. 

A desirable feature for such a digital simulation program is to have a built-in monitor of 
momenta and energy. The purpose of this appendix is to develop (in terms of previously 
identified state variables and system parameters) the expressions for total system angular 
and linear momentum vectors and the total system energy. 

The total angular momentum about the inertial reference can be expressed (from definition) 
as 


H = 



v) dm 


(D-1) 


with the summation over the number of bodies, NB, of the system, with x being the vector 
positioning the elemental mass, dm, from the inertial origin, with v being the absolute 
velocity of dm, and with integration taken over the volume of the body, Vj. Also, from 
definition, the total linear momentum with respect to the inertial frame is 

NB 

- E 

J=i 


Blank 


/ 


V dm 


(D-2) 


Pr6C6d\n9 
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Now, consistent with the notion of a body fixed-axis system and with a consistent velocity 
field assumed (paragraph II.B of the main text), it follows that, over the volume of the 
body. 




(D-3) 


and 


V = Vrj + w. X (p^ + r,)+ (D-4) 

On substituting equations D-3 and D-4 into D-1 and D-2 and integrating, it becomes clear 
that the first six elements of the product 

Id) ' ["’ll Ml 

’ (D-5) 

{Py \ 

fPJ f 

are projections of the body’s angular and linear momentum vectors onto the moving-body 
axis system. In fact, includes the effect of momentum wheels (equation 11-109 of the 
main text), which must be accounted for. Thus, the angular momentum of the j'* body 
about its body origin is 


H, = UJ {p„f, (D-6) 

whereas the linear momentum of the body is 

where le J is the unit-vector basis associated with the body fixed-reference triad. 



Rotation transformations that relate vector components in each body system to the inertial 
system exist, and position vector from the inertial origin to the reference point of each body 
exists. It follows that 






= E 

;=i 


loRj] iPvtj 


and that 


(D-8) 


E 

i= i 


( "j ♦ ’‘ri X 'l ) 



+ jsK*([oRj] Ipj)] 



(D-9) 


The total angular and linear momentum vectors are calculated by the program in the manner 
indicated in equations D-8 and D-9. For a variety of torque/force-free configurations that 
have been examined, momentum has been conserved within acceptable numerical tolerances. 

The total energy is calculated (equations 11-38 and 11-42 of the main text) as 


NB 

T + V = 

j=i 


([UJ. [ml3 \V\. + m. [k], ill.) 


(D-10) 


The kinetic energy contribution of embedded momentum wheels is included (as it must be) 
because [m]j includes momentum-wheel inertial coupling terms and jlJ|j includes momentum- 
wheel spin rates, (0^). 

Potential energy, in addition to that shown in equation D-10, comes about if there is a 
“sprung” hinge (e.g., associated with the coordinate). If the spring force/torque is linear 
with , additional potential energy is 

VcMBon.* = Si E (D-11) 

k 
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ROOT-LOCUS SOLUTION TECHNIQUES 


INTRODUCTION 

The root-locus solution procedure that is included in the digital program was supplied by 
Goddard Space Flight Center (GSFC) personnel and is included here for completeness. The 
following discussion has been extracted from GSFC Branch Report No. 254, dated October 2, 
1973 (Raymond V. Welch, author). 

Two expressions are identified as the basis for initiating the solution process; 


P(s) 
K— ^ 
Q(s) 


= 1 


(E-1) 



(E-2) 


Expressions E-1 and E-2 contain polynomial expressions for the conventional open-loop 
expression. 


KG(s) H(s) = KP(s)/Q(s) 


(E-3) 


The solution process requires a starting point (known to lie somewhere on the loci) from 
which to generate the desired locus; therefore, assume that a known value of s exists, (e.g., 
) so that 



(E-4) 


(i.e., Sjj is a point on the locus). A good starting point might be an open-loop pole or zero 
because these points are usually known a priori-, however, any point along the branch of the 
locus to be determined may be used. The locus is then traced, using the following procedure. 

E-3 
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Draw a “small” circle of radius, r, centered at s^^ , and define s^. to be the values of s which 
lie on this circle (figure E-1). These values of s^, are determined analytically by 


.9 

s = s + re^ 

C O 


(E-5] 


Segment Along One 
Branch of the Locus 


jw 




where 6 is measured as the angle generated by a counterclockwise rotation from the positive 
real axis of the s-plane. If the locus does not terminate inside the circle, at least two values 
of 6 must exist between zero and 360 degrees so that 



= 180° 


(E-6) 


(More than two values of d will exist that satisfy this equation if breakway points of the 
locus are encircled or if the circle is large enou^ to intersect other branches of the locus. 
Consequently, to avoid changing branches, r must be kept smaller than the distance between 
the branches of the locus.) If, for 0 = 0 j and 0 = $ 2 , the foregoing equation is satisfied, 

Sj.j and Sj .2 are points on the locus (figure E-1) where 


s , = 
cl o 


s 


c2 


s + re^®2 
O 


(E-7) 


These roots are found by an iterative process. 

STEP 2 

Define i//(0) by 



(E-8) 


where s^. is defined above and 0 is any arbitrary angle. If i//(0) does not equal 180°, increment 
0 by A0 and reevaluate i|/(0). Continue this process until i//(0) passes through 180°. Because 
i//(0) is a monotonic function across the locus, i//(0) crossing 180° implies that the locus has 
been crossed. When i|/(0) passes through 180°, redefine A0 as 


AO, 


= -KA0 


old 


0<K<1 


(E-9) 


and again calculate \l/(0). If ^1^(0) does not equal 1 80° for this point, increment 0 by this new 
A0 and recalculate i//(0). Continue this process until i//(0) again crosses 180°. Reduce A0 
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again by the foregoing process, and repeat the previous operation until Ad = 6, where 6 is 
some predefined small positive number. At this time, \}j{d) will be 180 ±e°, where e is a 
“small” angle whose size is a function of 5 ; thus, a root on the locus has been found that 
is either s^,^ or Sg 2 ) depending on the initial choice of 8 and on the initial sign of Ad. For 
this subprogram, the initial value of 0 is optional, but, if no choice is made, a value of 1 80° 
is used. Also, the initial increment. Ad, is set at 10° with the change in Ad for every 180° 
crossing defined by 


A0 


new 


1 

To 


AO 


old 


(E-10) 


It was found that 5 = 10"^ was small enough to ensure e < 0.001° except for possibly 
values of s near the open-loop poles or zeros. Assume that d is chosen initially so that the 
root found is located at the point, s^j , as previously defined and as shown in figure E-1, and 
suppose that it is desired to continue the locus in this direction. 

STEP 3 

Draw a circle of radius, £ , centered at s^j . Again, if the locus does not terminate inside the 
circle, the locus will intersect the circle in at least two points. Because the circles located 
at Sq and s^j are of equal radius, one of these points is s^ as shown in figure E-2. This root 
is eliminated from the search routine by restricting the range of 8 to 

- 120° < 0 < + 120° (E-11) 

These limits are chosen because they are the points at which the circles located at s^^ and 
Sj,j intersect. Thus, the search for a new root is conducted only in a previously unsearched 
region. The initial angle is chosen at 0 = - 120° A0>o, and Step 2 is repeated until 

A0 = S (i.e., another root is found). 

STEP 4 

Draw a circle of radius, r, centered at the root found in the previous step, and repeat Step 2 
with the restrictions on 0 as defined in Step 3. 

Repeat Step 4 to determine the roots along one branch of the locus. The value of the gain, 
K, for each of these roots is calculated from 


K = 


Q(s) 

P(s) 


(E-1 2) 
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a 

Figure E-2. Search area definition for finding roots 
after the first root is found. 
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TYPICAL KULTIBODY CONFIGURATION 



Note 

1 . 


All oomentun wheels must 
also have « sensor point, 

Body t must always be 
positioned relative to 
the Inertial reference. 
Hinge 1 is always between 
Body I and the inertial 
reference. 


Lesend: 

9 

Body Reference 

Point 

• 

Hinge Point 


0 

Sensor Point 



Komeotum Wheel 



Inertial Reference 



IK . 

1 





TYPICAL TWO BODY/HIHGE POINT ARRANGSHENT 



Preceding Page Blank 
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PROGRAM SIMULATION INTEGER ABRAYS 



Body No 


Topology, Size • 2 x NH 
1 2 3 4 5 

Body n has frame q 
Body m has frame p 


Example: 



Refer to Column 3 (Hinge 3) and 
note Body 4 Interfaces with Body 
1 via Hinge 3. 

No. of Elastic Hodes/Body. Size ■ 1 x NB 


Example: Body No. 4 has 2 elastic modes, 

Body No. 2 has no elastic modes. 


I IFTSMW I Sensor Point Locations, Size » 1 x NS 

Sensor 2 3 4 5 6 

[3 4 2 4 4 3l 


IHDATA 


Exan^le: Sensor Point is located on Body 

(as are Sensor Points and 

Hinge Information, Size * 7 x NH 

1 2 3 4 5 Hinge No. 


mpE 

2 1 2 2 U 


Note: 

h 

0 0 0 1 0 


1, ITYPE - Euler Type 

02 

0 0 0 1 2 


2. Elements in rows 2-7 identify 

03 

0 0 0 0 1 


hinge constraint type. 

x 

0 1111 


0: Free/forced 

. 



1: Fixed Constraint 

y 

0 10 11 


(Zero relative velocity) 


0 1110 


2: Rheonomlc Constraint 


1 


(time dependent) 


Additional Remarks ; 

1. The number of Betas in the state vector 
equals the sum of "zeros" plus the sum 
of the "twos" (excluding row 1} in the 
array. 

2. The number of Lambdas (constraints) 
equals the number of "nonzeros" 
(excluding row 1) in the array. 



12 3 


AMO 


5 14 

2 13 

110 


0] 02 03 
Jl Jz J} 


Momentum Wheel Information, Size » 3 x NOFMO and 

2 X NOFMO 

Homentum Wheel No. 

Sensor Point No. 

Spin Axis.(l, 2, or 3) 

1 = Active; 0 » Constant Speed 
Wheel Rates (Initial) 

Wheel Inertia about Spin Axis 






EULER AMCLE PEtaOTATIOS CAMDIDATES i. ITYPE 




ITYPE - 


1 

2 

3 

4 

S 

6 

7 


9 

10 

11 

12 

1st 

Kotatlon 

Axis of 


1 

1 

1 

T 

rr 

2 

2 


3 

3 

3 

3 

2nd 

Rotation 

Axis of 

92 

2 

2 

3 



3 

1 


1 

1 

2 

2 

3rd 

Rotation 

Axis of 

Bj 

3 

1 

1 


li 

2 

2 

3 

2 

3 

3 

1 


Exajcple: Using ITTfPE - 6 



CORSOLIDATION OF KINEHATICAL COEFFICIENTS 
(The "b" Coefficients) 




fU)3 

{U)^ 


(B), 

( 8)2 

< !> 

(S>3 

(3)c 


Note : 

1. Only a single entry Is in partition for 
Hinge 1 and it is a q type parti-* 

tlon froa Che connection klnematlcal 
array above, 

2. All other row partitions have both a p 

and a q constituent for each "hinge" | 
between contiguous bodies. I 
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CONNECTION KINEMATICS ^ TYPICAL HINGE 


a "p" Point 


a *'q** Point 




A62 


A63 

< . 


fixi 


Ax2 


Ax3 



-Tt“l R R 
q p p m 


R R 0 
q p p m p 

R 

q n 


' 

R 0 
q n q 

- R S 
p ffl mp 

- R 
p m 

- R h 
p m p 

R R S 
p q q n nq 

R R 
p q q n 

R R h 
p q q n q 


!• Is a rotation transformation from 
^ the J to the 1 reference frame. 

2. In general, 

is a transformation relating Euler 
rates to body axes projections of 
the angular velocity vector, 

is a modal slope matrix in the i 
frame , 

h. is a modal deflection matrix in the 
1 frame , 

S. is a skew symmetric matrix whose 
^ elements are the components of the 
vector 1 to J in the ith frame. 


Pil- 


V 


l”n) 

j?) I 

I ml 


STATE VECTOR ARRANGEMENT 


CONSTRAINT FORCE/TORQUE VECTOR 


U) - 


^1 

^2 

^5 

^7_ 

la 

AlO 

^11 

^12 

^13 


'V3 

(Ti>4 

(T2)^ 

(T2)s 

(T3>5 

<Vs 


(F refers to Hinge i 


ty) 


'{u>r 


'{U)i 

{U)2 


{U>2 

{U >3 


{U )3 

<u)i. 

; (y) - 

{6)4 

u’h 


'u’h 

U)2 


U)2 

U )3 


U >3 

{«>4 


{C >4 

’ai’ 


si" 

82 


82 

83 


63 

84 


84 

85 


65 

86 


06 

67 

► < 

87 

8g 


88 

89 


89 

810 


8)0 

811 


011 

812 


012 

813 


013 

014 


8)4 

8)5 


01$ 

8|6 


016 

837 


-11- 

'4l" 


ly 

82 


82 

«3 


«3 

N > 




where (U) 




-For this example, prescribed 
as function of time, l.e., a 
rheonomlc constraint as 
noted In IHDATA 
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